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Preface 



Much interest has been drawn in recent years to the concept and realiza- 
tion of Nanoelectromechanical systems (NEMS). NEMS are nanoscale devices 
that combine mechanical and electrical dynamics in a strong interplay. The 
shuttle devices are a particular kind of Nanoelectromechanical systems. The 
characteristic component that gives the name to these devices is an oscillating 
quantum dot of nanometer size that transfers electrons one-by-one between 
a source and a drain lead. The device represents the nano-scale analog of an 
electromechanical bell in which a metallic ball placed between the plates of 
a capacitor starts to oscillate when a high voltage is applied to the plates. 

This thesis contains the description and analysis of the dynamics of two 
realizations of quantum shuttle devices. We describe the dynamics using the 
Generalized Master Equation approach: a well-suited method to treat this 
kind of open quantum systems. We also classify the operating modes in three 
different regimes: the tunneling, the shuttling and the coexistence regime. 
The characterization of these regimes is given in terms of three investigation 
tools: Wigner distribution functions, current and current-noise. The essen- 
tial dynamics of these regimes is captured by three simplified models whose 
derivation from the full description is possible due to the time scale separa- 
tion of the particular regime. We also obtain from these simplified models a 
more intuitive picture of the variety of different dynamics exhibited by the 
shuttle devices. 
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Chapter 1 
Introduction 



In this chapter we give a short introduction to the world of nanoelectrome- 
chanical systems. We then focus our attention on a particular kind of device 
called electron shuttle. We sketch the basic operating regime and give an 
overview of the different theoretical models that have been proposed to de- 
scribe the dynamics of such devices. We report on the two main realizations 
of shuttle devices and close the chapter with an outline of the contents of 
this thesis. 



1.1 NEMS 

Much interest has been drawn in recent years on the concepts and realiza- 
tion of Nanoelectromechanical systems (NEMS). NEMS are nanoscale devices 
that combine mechanical and electrical dynamics in a strong interplay. This 
property makes them interesting both from a technological and fundamental 
point of view. They are extremely sensitive mass and position detectors. 
Due to their very high mechanical frequency one can even think of using 
them as the basis for new form of mechanical computers. From the point of 
view of fundamental research they represent extremely good tools to probe 
directly the basic quantum mechanical laws. They could represent the first 
man-made structures on which the mechanical zero point fluctuation can be 
detected. They also rise the question on the limiting dimension for persis- 
tence of mechanical coherence. In general one of the fascinating aspects of 
these objects is their mesoscopic character: they share with the macroscopic 
world the large number of atoms of which they are made (typically of millions 
of atoms) but on the other hand their behaviour is (or should be) already 
significantly determined by quantum mechanics. 
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1.2 A new transport regime 

The shuttle devices are a particular kind of NEMS. The characteristic compo- 
nent that gives the name to these devices is an oscillating object of nanometer 
size that transfers electrons one-by-one between a source and a drain lead. 
The device represents the nano-scale analog of an electromechanical bell in 
which a metallic ball placed between the plates of a capacitor starts to os- 
cillate when a high voltage is applied to the plates. The oscillations are 
sustained by the external bias that pumps energy into the mechanical sys- 
tem: when the ball is in contact with the negatively biased plate it gets 
charged and the electrostatic field drives it towards the other capacitor plate 
where the ball releases the electrons and returns back due to the oscillator 
restoring forces 1 and the cycle starts again. 

In the first proposal [1] of a shuttle device the movable carrier is a metallic 
grain confined into a harmonic potential by elastically deformable organic 
molecular links attached to the leads. The transfer of charge is governed by 
tunneling events, the tunneling amplitude being modulated by the position of 
the oscillating grain. The exponential dependence of the tunneling amplitude 
of the grain position leads to an alternating opening and closing of the left 
and right tunneling channels that resembles the charging and discharging 
dynamics of the macroscopic analog. 

Different models for shuttle devices have been proposed in the literature 
since this first seminal work by Gorelik et al. [1]. The mechanical degree 
of freedom has been treated classically (using harmonic [2, 3, 4, 5] or more 
realistic potentials [6]) and quantum mechanically [7, 8, 9]. Armour and 
MacKinnon proposed a model with the oscillating grain flanked by two static 
quantum dots [7, 10, 11, 12]. More generally the shuttling mechanism has 
been applied to Cooper pair transport [13, 14] and pumping of superconduct- 
ing phase [15] or magnetic polarization [16]. 

The essential feature of the nano-scale realization is the quantity trans- 
ferred per cycle (a charge up to 10 10 electrons for a macroscopic bell) that 
is scaled down to 1 quantum unit (electron, spin, Cooper pair in the differ- 
ent realizations). We can already guess the basic properties of the shuttle 
transport: 

1. Charge-position correlation: the shuttling dot loads the charge on one 
side and transfers it on the other side, it releases it and returns back 

1 Due to the large amount of electrons in this macroscopic realization the ball gets 
positively charged at the second plate by loosing some extra electrons and the restoring 
force contains also an electrostatic component. The system is perfectly symmetric under 
commutation of charge sign. 
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to the starting point; 

2. Matching between electronic and mechanical characteristic times (non- 
adiabaticity); 

3. Quantized current determined by the mechanical frequency; 

4. Low current fluctuations: the stochasticity of the tunneling event is 
suppressed due to an interplay between mechanical and electrical prop- 
erties. The tunneling event is only probable at some particular short 
time periods fixed by the mechanical dynamics (i.e. when the oscillat- 
ing dot is close to a specific lead). 

1.3 Experimental implementations 

An experimental realization of the shuttle device has been produced by Erbe 
et al. [17]. The structure consists of a cantilever with a quantum island at 
the top placed between source and drain leads. Two lateral gates can set 
the cantilever into motion via a capacitive coupling. An ac voltage applied 
at these gates makes the cantilever vibrating and brings the tip alternatively 
closer to the source or drain lead and thus allows the shuttling of electrons. 

The device (shown in Fig. 1.1) is built out of silicon-on-insulator (SOI) 
materials (using Au for the metallic parts) using etch mask techniques and 
optical and electron beam lithography. The cantilever is lfim long and the 
corresponding resonant frequency is of the order of 100 MHz. The source elec- 
trode and the cantilever are at an average distance of approximately 0.1 \im. 
Shuttling experiments have been performed by Erbe et al. at different tem- 
peratures. For experiments at 4.2 K and 12 K they measured a pronounced 
peak in the current through the cantilever for a driving frequency of approx- 
imately 120 MHz corresponding to the natural frequency of the first mode 
of the oscillator. The peak in the current corresponds to a rate of shuttle 
success of about one electron per 9 mechanical cycles. The Erbe experiment 
is very close to the original proposal by Gorelik. The only difference is in the 
external driving of the mechanical oscillations. In the original proposal the 
bias was time independent and the driving induced by the electrostatic force 
on the charged oscillating island. The initially stochastic tunneling events 
would eventually cause the shuttling instability and drive the system into 
a self-sustained mechanical oscillation combined with periodic charging and 
discharging events. 

Another experiment often mentioned in the context of quantum shuttles 
is the Ceo-experiment by Park et al. [18] In this experiment a Ceo molecule is 
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Figure 1.1: Electron micrograph of the nano-mechanical resonator. The cantilever 
(C) can be set into motion by applying an ac-voltage between the two gates Gl and 
G2, and by applying a bias across the metallic tip of the cantilever (the island) 
electrons are shuttled from the source (S) to the drain (D). The picture is taken 
from [17]. 




Distance 



Figure 1.2: The Cqq experiment by Park et al. The Cqq molecule (of mass M) can 
be considered as being attached to a spring with spring constant k and corresponding 
quantized excitation energy hf of the order of 5 meV. The figure is taken from [18]. 
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deposited in the gap between two gold electrodes. The gap, produced with 
break junction technique, has a width of 1 nm. The molecule (of diameter 0.7 
nm) is bound to the electrodes due to van der Waals interaction. Around the 
equilibrium position the potential can be approximated by a harmonic po- 
tential and the molecule can be considered as attached by springs. We show 
a schematic representation of this idea in Fig. 1.2. In the experiment Park 
et al. sweep the voltage across the junction and register sudden increases in 
the current. The steps are separated by 5 meV. Since the lowest internal ex- 
citation energy of the Cqq molecule is 35 meV one concludes that the slower 
center of mass motion could be involved in the process. This hypothesis is 
confirmed by the fact that the separation between the steps is independent 
of the charge on the C 60 molecule and the theoretical estimate for the energy 
corresponding to the center of mass oscillation in the van der Waals potential 
is exactly 5 meV. The IV-curve measured in these experiments can be inter- 
preted in terms of shuttling [19], but also alternative explanations have been 
promoted [20, 21, 22]. Whether we are in presence of coherent or incoherent 
shuttling transport and to what extent the shuttling mechanism is involved 
in this set-up cannot be completely clarified by current measurement. The 
current-current correlation (also called current noise) was proposed by Fe- 
dorets et al. [19] for a better understanding of the underlying dynamics of 
the current jumps. For this reason the calculation of noise in shuttle devices 
has been performed by many different groups. We proposed a completely 
quantum formulation [23] that is also explained in detail in this thesis. 

1.4 This thesis 

This thesis contains the description and analysis of the dynamics of two 
realizations of quantum shuttle devices. The models we consider describe 
both the mechanical and electrical degrees of freedom quantum mechanically. 
For the single dot quantum shuttle we extended an existing classical model 
proposed by Gorelik et al. [1]. For the triple dot case we adopted the model 
invented by Armour and MacKinnon [7]. In the following we outline the 
contents of the thesis: 

In Chapter 2 we introduce the two models called Single Dot Quantum 
Shuttle and Triple Dot Quantum Shuttle, the first being the quantum ex- 
tension also for the mechanical degree of freedom of the model originally 
proposed by Gorelik et al. [1] while the second is the model invented by 
Armour and MacKinnon [7]. Also in this model the oscillator is treated 
quantum mechanically and the central moving dot is flanked by two static 
dots. 
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We dedicate Chapter 3 to the derivation of the Generalized Master Equa- 
tion (GME) that describes the shuttle device dynamics. Due to the different 
coupling strengths we treat the mechanical and electrical baths with two dif- 
ferent approaches. The Gurvitz approach for the electrical and the standard 
Born-Markov approximation for the mechanical bath. The derivation a la 
Gurvitz represents a large part of this chapter. In order to facilitate the 
understanding of this non-standard method and appreciate our extension to 
the shuttle device we have given a long introduction in which we analyze in 
great detail simpler models with increasing physical complexity. This shows 
on one hand the essence of the derivation in simpler cases and also under- 
lines the potentiality of this approach. An important aspect of this method 
is also to be a natural prelude to full counting statistics since it naturally 
produces a GME that counts the number of electron which tunneled through 
the device at a certain time. We close the chapter with the description of 
the numerical iterative method that we adopted for the calculation of the 
stationary solution of the GME. 

In Chapter 4 we introduce the concept of Wigner distribution function 
and derive the corresponding Klein-Kramers equation for the SDQS start- 
ing from the GME that we obtained in the previous chapter. The Wigner 
function description is motivated by the effort to keep the complete quantum 
treatment we achieved with the GME without losing as much as possible the 
intuitive classical picture and with the possibility to handle the quantum- 
classical correspondence. 

Chapter 5 is dedicated to the definition and application of the three inves- 
tigation tools we have chosen to analyze the properties of shuttle devices: the 
charge resolved phase-space distribution, the current and the current-noise. 
The phase space analysis reveals the shuttling transition and the charge po- 
sition correlation typical of this operating regime. It also gives a very clean 
way to appreciate "geometrically" the quantum to classical transition of the 
shuttling behaviour for different device realizations. The second investigation 
tool that we consider is the current. From the current calculation we obtain 
also in the quantum treatment the quantized value of one electron per cycle 
found in the semiclassical treatments of similar devices. We then present 
current-noise calculations based on the MacDonald formula. The derivation 
is strongly dependent on the derivation a la Gurvitz of the n-resolved GME. 
The low noise quasi-deterministic behaviour of the shuttling transport is clear 
from the extremely low Fano factors found for this regime. In general we are 
able with all the three investigation tools to identify three operating regimes 
of shuttle devices: the tunneling, shuttling and coexistence regimes. 

Chapter 6 is dedicated to a qualitative description of the these regimes 
and also to the identification of the relation between different times and 
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length scales that define the three regimes in terms of the model parameters. 

In Chapter 7 we consider separately the tunneling, shuttling and coexis- 
tence regime. The specific separation of time scales allows us to identify the 
relevant variables and describe each regime by a specific simplified model. 
Models for the tunneling, shuttling and coexistence regime are analyzed in 
this chapter. We also give a comparison with the full description in terms of 
Wigner distributions, current and current-noise to prove that the models, at 
least in the limits set by the chosen investigation tools, capture the relevant 
dynamics. 

A summary of the arguments treated in this thesis opens Chapter 8. We 
conclude with a list of some of the open questions that could encourage a 
continuation of the present work. 
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We describe in this chapter two models of quantum shuttles: the Single- 
Dot Quantum Shuttle and the Triple-Dot Quantum Shuttle. Due to the 
nanometer size of these devices we decide to treat quantum mechanically 
not only the electrical but also the mechanical dynamics. This approach was 
suggested by the work of Armour and MacKinnon [7] for the triple dot device 
and implemented for the first time by us in the single dot device. 

2.1 Single-Dot Quantum Shuttle 

The Single-Dot Quantum Shuttle (SDQS) consists of a movable quantum 
dot (QD) suspended between source and drain leads. One can imagine the 
dot attached to the tip of a cantilever or connected to the leads by some 
soft legands or embedded into an elastic matrix. In the model the center of 
mass of the nanoparticle is confined to a potential that, at least for small 
diplacements from its equilibrium position, can be considered harmonic. We 
give a schematic visualization of the device in figure 2.1. 

Due to its small diameter, the QD has a very small capacitance and thus 
a charging energy that exceeds (in the most recent realizations almost at 
room temperature [24]) the thermal energy ksT 1 . For this reason we assume 
that only one excess electron can occupy the device (Coulomb blockade) 
and we describe the electronic state of the central dot as a two-level system 
(empty/charged). Electrons can tunnel between leads and dot with tunneling 
amplitudes which are exponentially dependent on the position of the central 
island. This is due to the exponentially decreasing/increasing overlapping of 

*A quick estimate of the charging energy can be obtained for an isolated 2D metallic 
disk: e 2 /C = e 2 / (8e r eoR) where R is the disk radius and e r = 13 in GaAs. For a dot of 
radius 10 nm this yields e 2 /C = 20meV = /c B 230K 
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1-10 nm 



Figure 2.1: Schematic representation of the Single-dot Shuttle: electrons tunnel 
from the left lead at chemical potential (hl) to the quantum dot and eventually to 
the right lead at lower chemical potential [ir. The position dependent tunneling 
amplitudes are indicated. X is the displacement from the equilibrium position. The 
springs represent the harmonic potential in which the central dot can move. 



the electronic wave functions. 

The Hamiltonian of the model reads: 

H = H sys + H + -f/bath + H tun + H mt (2.1) 

where 

H svs = 1 — muj 2 x 2 + (si — e£x)c\ci 

y 2m 2 V 711 

Pleads = J2( £l * C l C h + £ rA k Cr k ) 

k (2.2) 

#tun = y^jT^c}^ + T r {x)cl k c x ] + h.c. 

k 

Hbath + H [nt = generic heat bath 

Using the language of quantum optics we call the movable grain alone the 
system. This is then coupled to two electric baths (the leads) and a generic 
heat bath. The system is described by a single electronic level of energy 
E\ and a harmonic oscillator of mass m and frequency u. When the dot is 
charged the electrostatic force (e£) acts on the grain and gives the electrical 
influence on the mechanical dynamics. The electric field S is generated by 
the voltage drop between left and right lead. In our model, though, it is 
kept as an external parameter, also in view of the fact that we will always 
assume the potential drop to be much larger than any other energy scale 
of the system (with the only exception of the charging energy of the dot). 
The operator form x,p for the mechanical variables is due to the quantum 
treatment of the harmonic oscillator. In terms of creation and annihilation 
operators for oscillator excitations we would write: 



20 



2.1. SINGLE-DOT QUANTUM SHUTTLE 




The leads are Fermi seas kept at two different chemical potentials (//£, 
and nn) by the external applied voltage (AV = — fiR)/e ). The oscillator 
is immersed into a dissipative environment that we model as a collection of 
bosons and is coupled to that by a weak bilinear interaction: 

#bath = ^ hWqdq^dq 

(2.4) 

q 

where the bosons have been labelled by their wave number q. The damping 
rate is given by: 

7(w) = 2ng 2 D(iu) (2.5) 

where D(u) is the density of states for the bosonic bath at the frequency of 
the system oscillator. A bath that generates a frequency independent 7 is 
called Ohmic. 

The coupling to the electric baths is introduced by the tunneling Hamilto- 
nian H tmi . The tunneling amplitudes Ti(x) and T r (x) depend exponentially 
on the position operator x and represent the mechanical feedback on the 
electrical dynamics: 



Ti t r(x) = ti, r exp(=F£/A) (2.6) 

where A is the tunneling length. The tunneling rates from and to the leads 
(f l,r) can be expressed in terms of the amplitudes: 

r L , R = (r L , R m = (y Dl ' rGX v ( t t) '^ |2 ) (2 - 7) 

where D L R are the densities of states of the left and right lead respectively 
and the average is taken with respect to the quantum state of the oscillator. 

The model presents three relevant time scales: the period of the oscillator 
2n/u, the inverse of the damping rate I/7 and the average injection/ejection 
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time l/V L>R . It is possible also to identify three important length scales: 

the zero point uncertainty Ax z = \J the tunneling length A and the 

displaced oscillator equilibrium position d = Different relations between 
time and length scales distinguish different operating regimes of the SDQS. 



2.2 Triple-Dot Quantum Shuttle 

The Triple-Dot Quantum Shuttle (TDQS) was proposed by Armour and 
MacKinnon [7]. The system consists of an array of three QD's: a movable 
dot, that we assume confined to a harmonic potential, flanked by two static 
ones. Relying on low temperature and on the low capacitance of the system 
with respect to the leads, we again assume strong Coulomb blockade: only 
one electron at a time can occupy the three-dot device. The Hamiltonian for 
the model reads: 

H = H sys + Hi cads + //bath + H tun + H mt (2.8) 

Only the system and tunneling part of the Hamiltonian differ from the 
one dot case: 

H sys =e |0)(0| + ^f\L)(L\ - ^\C)(C\ - ^f\R)(R\ + tiu; (dU + 

+ t R (x)(\C)(R\ + \R)(C\) + t L (x)(\C)(L\ + \L)(C\) 

//tun =5>( C ;jO)<L| +c lk \L)(0\) +T r (4 k \0)(R\ + c rk \R){0\) 
k 

(2.9) 

where H osc is the harmonic-oscillator Hamiltonian, \a), a = 0,L,C,R are the 
vectors that span the electronic part of the system Hilbert space . The tun- 
able injection and ejection energies (the energy levels of the outer dots, that 
we can assume fixed by external gates) simulate a controlled bias through the 
device (AV) and the position dependent tunneling amplitudes are now be- 
tween elements within the system. These amplitude are assumed to be expo- 
nentially dependent on the position of the central dot ti(x) = —V G e~( Xa+x ^ x 
and tn(x) = — Voe - ^ ^^. Tunneling from the leads is allowed only to the 
nearest dot and the corresponding tunneling amplitude is independent of the 
position of the oscillator. The "device bias" AV also gives rise to an electro- 
static force on the central dot, when charged. A schematic representation of 
the Triple-dot Shuttle is given in figure 2.2. 

For reasons that will become clearer in the following, we assume that all 
the energy levels of the system (except the Coulomb charging energy that 
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1-10 nm 



Figure 2.2: Schematic representation of the Triple-dot Shuttle: the leads and the 
three-dot array are represented. The arrows mimic the electrical dynamics. Single 
and double arrows indicate that the tunneling from and to the lead is always in a 
given direction and incoherent while the internal dynamics of the system is subject 
to coherent oscillations. The mechanical motion of the central dot confined to a 
harmonic potential is represented by the springs. 



ensures the strong Coulomb blockade regime) lie well inside the bias window. 
In practice we will take the limit //£ — > oo and jiR — > — oo. This is reflected 
in the directional flow of electrons from the source and to the drain. 
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Chapter 3 

Generalized Master Equation 



The state of a physical system is determined by the measurement of a certain 
number of observables. Repeated measurements of a given observable always 
return the same expectation value when the system is in an eigenstate for that 
particular observable. The uncertainty principle ensures us that, for quantum 
systems, there are incompatible observables that can not be measured at the 
same time with indefinite precision. 

Given a generic quantum system S and a complete set of compatible ob- 
servables Ai [25], an eigenstate of the system for all observables is defined by 
the set of the corresponding expectation values, i.e. the quantum numbers 
a,{. Each of the possible sets of expectation values is associated with an eigen- 
vector in the Hilbert space of the system. More precisely the Hilbert space 
of the system is spanned by the eigenvectors of a complete set of compatible 
observables. 

A pure state of the system is represented by a radius (class of equivalence 
of normalized vectors with arbitrary phase) of this Hilbert space. We call 
a representative vector of a radius. Observables are associated to Hermitian 
operators on the Hilbert space of the system. The dynamics of the quantum 
system is governed by the Hamiltonian operator, i.e. the operator associated 
to the observable energy. Given an initial vector, the Schrodinger equation 
prescribes the evolution of this vector at all times: 

ihj t m)) = n\m) (3.i) 

with the initial condition | -0(0) ) = \ipo). An equivalent formulation of the 
dynamics can be given in terms of projector operators \ip){if}\. A projector is 
independent from the arbitrary phase of the vector it is then equivalent 
to a radius of the Hilbert space and represents a pure state of the system. 
Using the Leibnitz theorem for derivatives and the Schrodinger equation we 
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derive the equation of Liouville-von Neumann: 



dp 
~dt 




(3.2) 



where p = \ip){if}\, [A, B] = AB — BA is the commutator of the operators A 
and B. The operator p is usually called density operator. For each basis of 
the Hilbert space all the operators have a matrix representation. The matrix 
that corresponds to the density operator is called density matrix. Each vector 
of the basis of the Hilbert space corresponds to a particular eigenstate of 
the system defined by a set of quantum numbers. The diagonal elements 
of the density matrix are called populations. Each population represents 
the probability that the system in the pure state I^X^I is in the eigenstate 
defined by the corresponding set of quantum numbers. The trace of the 
density matrix is one and supports this probabilistic interpretation. The 
off-diagonal terms of the density matrix are the coherencies of the system. 
They reflect the linear structure of the Hilbert space. A linear combination 
of eigenvectors gives rise to a pure state with non-zero coherencies. 

Not all density operators correspond to pure states. A convex linear 
combination of pure states iV'nXV'nl, n — lj ■■■,N is called statistical mixture: 



where P n G [0,1), ^2 n P n — 1- This is an incoherent superposition of pure 
states. Also statistical mixtures obey the Liuoville-von Neumann equation 
of motion (3.2). 

The master equation is an equation of motion for the populations. It is 
a coarse grained 1 equation that neglects coherencies. It was derived the first 
time by Pauli under the assumption that coherencies have random phases in 
time due to fast molecular dynamics. It reads: 



where P n is the population 2 of the eigenstate n and T nm is the rate of prob- 
ability flow from eigenstate n to m [26]. 

1 In the sense that it describes the effective dynamics on a time scaie long compared to 
the typical times of the fastest processes in the physical system. 

2 Since a density matrix without coherencies is a statistical mixture of eigenstates we 
have adopted the notation P n = p nn 



N 




(3.3) 



n=l 




(3.4) 



m 
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3.1 Coherent dynamics of small open systems 

The master equation is usually derived for models in which a "small" system 
with few degrees of freedom is in interaction with a "large" bath with effec- 
tively an infinite number of degrees of freedom. The Liouville von-Neumann 
equation of motion for the total density matrix is very complicated to solve 
and actually contains too much information since it also takes into account 
coherencies of the bath. It is useful to average it over bath variables and 
obtain an equation of motion for the density matrix of the system (the re- 
duced density matrix). With no further simplification this equation is called 
Generalized Master Equation (GME) since it involves not only the pop- 
ulations but also the coherencies of the small subsystem. The derivation 
of the GME from the equation of Liouville-von Neumann is far from trivial 
and also non- universal: it involves a series of approximations justified by 
the physical properties of the model at hand. Despite the apparent similar- 
ities, the two equations are deeply different: the equation of Liouville-von 
Neumann describes the reversible dynamics of a closed system; the GME, 
instead, describes the irreversible dynamics of an open system that continu- 
ously exchange energy with the bath 3 . 

Shuttle devices are small systems coupled to different baths (leads, ther- 
mal bath) but they maintain a high degree of correlation between electrical 
and mechanical degrees of freedom captured by the coherencies of the reduced 
density matrix. The GME seems to be a good candidate for the description 
of their dynamics. 

In the next two sections we will derive two GMEs using two different 
approaches. They are both necessary for the description of the shuttling 
devices since they correspond to the different coupling of the system to the 
mechanical and electrical baths. 

3.2 Quantum optical derivation 

The harmonic oscillator weakly coupled to a bosonic bath is a typical problem 
analyzed in quantum optics. This model well describes in shuttling devices 
the interaction of the mechanical degree of freedom of the NEMS with its 
environment. Following section 5.1 of the book "Quantum Noise" by C. W. 
Gardiner and P. Zoller [27] we start considering a small system S coupled to 
a large bath B described by the generic Hamiltonian: 

3 How can irreversibility be derived from reversibility? The solution of this dilemma lies 
in time scales: system+bath recurrence time is "infinite" on the time scale of the system. 
The GME holds on the time scales of the system. 
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H = Hs + H B + Hi 



(3.5) 



where H$ and H B respectively describe the dynamics of the decoupled system 
and bath and Hi represents the interaction between the two that we assume 
weak. The density operator p (the state of the system+bath) satisfies, in the 
Schrodinger picture, the equation of Liouville-von Neumann: 



p=-j_[H s + H B + H h p\ (3.6) 
The state of the system is described by the reduced density matrix a: 



a = Tr B {p} 



(3.7) 



where Tr B indicates the partial trace over the bath degrees of freedom. Our 
task is to derive from (3.6) an equation of motion for the reduced density 
matrix a. 



3.2.1 Interaction picture 

We start by going to the interaction picture and we use as non-interacting 
Hamiltonian H S + H B . We indicate all the operators in the interaction picture 
with a tilde. The total density operator in the interaction picture reads: 



p(t) = exp -{H s + H B )t 
In 

and obeys the equation of motion: 



p(t) exp 



--(H s + H B )t 

h 



(3.8) 



m = -\mt),p{t)} 



where 



(3.9) 



Hit) = exp 



h 



(H s + H B )t 



Hi exp 



--(H s + H B )t 



From (3.7) and (3.8) it follows that 



(3.10) 



a{t) = Tr B <^ exp 



h 



[H s + H B )t 



pit) exp 



h 



[H s + H B )t 



(3.11) 



The exponentials of H B can be cancelled using the cyclic property of the 
trace since H B depends only on bath variables. We get 
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a(t) = exp 



--H s t 
n 



a {t) exp 



n 



(3.12) 



where 



a{t) = Tr B {p(t)} 



(3.13) 



In other terms the interaction picture for the reduced density matrix is ef- 
fectively obtained only from the non-interacting Hamiltonian for the system 



3.2.2 Initial conditions 

We assume that the system and the bath are initially independent, the initial 
total density operator p is then factorized into the tensor product: 



p(0) = <t(0) ®p B 
For definiteness we assume the bath to be in thermal equilibrium: 



PB 



Tree-^B 

where j3 = 1/ksT is the inverse temperature. 



(3.14) 



(3.15) 



3.2.3 Reformulation of the equation of motion 

It is most important for the derivation of the GME to recast the original 
equation of motion in the interaction picture (3.9) into an integro-differential 
form. The integral from to £ of (3.9) reads 



and inserted back into (3.9) itself gives 



Pit) = -![£:(*), p(0)] dt'H.it), mt'),p(t') 



(3.16) 



(3.17) 
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3.2.4 Average over the bath variables 

We take the partial trace over the bath variables on both sides of (3.17). 
From the definition of the reduced density operator (3.7) we obtain, for the 
LHS, a. We assume that the first term in the RHS vanishes, namely 

Tr B [#i,p(0)] = (3.18) 

where p(0) = p(0) = o"(0)(g>pB- This means that the interaction Hamiltonian 
has a bath component with zero average. It is not difficult to fulfill this con- 
dition in general by a redefinition of the system and interaction Hamiltonian 
that subtracts the average of the bath component from the latter. 

3.2.5 Weak coupling 

We assume that H\ is only a small perturbation of H$ and H B . This condition 
allows a factorization at all times of the total density operator into its system 
and bath components. The density operator of the bath is also taken as 
constant in time: 

p(t)^a(t)®p B (3.19) 

The factorization assumption can be weakened. We introduce for this 
purpose the notion of correlation function. Given a physical system in the 
state described by the stationary density operator p and two operators 0\ 
and 2 the correlation function between 0\ and 2 in this order and at times 
t and t' is: 

C 0l o 2 (M') = T^pOiWC^O} (3.20) 

where the trace is taken over the the Hilbert space of the system and the 
operators are in Heisenberg picture. Returning to our system-bath model, 
we assume that the interaction Hamiltonian is a sum of operators in the form 
-Fs^b- The minimal requirement for the weak coupling approximation is that 
the correlation functions of the bath are not influenced by the state of the 
system. In formulas: 

Tr B {[A B (t), [A B (t'),p(t')]]} « a{t') ® Tr B {[A B (t), [A B (t'), p B ]]} (3.21) 

3.2.6 Markov approximation 

The integro-differential equation for a obtained in the weak coupling approx- 
imation is non-local in time. The state of the system at time t depends on 
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the history of the model starting from the initial time t — 0. This is the 
meaning of the integral on the RHS of the equation 4 

1 I* ~ 

* = ^'Tr B {[#i(0, [#i(*V(0 ®pb]]}. (3.22) 

obtained from (3.21) by tracing over bath variables. Due to the different 
sizes and the weak coupling the effects on the bath of the interaction with the 
system are negligible. The bath is a classical macroscopic object in thermal 
equilibrium. Its stationary state is a thermal state: an incoherent statistical 
mixture of energy eigenstates. The coherencies in the bath introduced by the 
interaction with the system decay on a time scale called the correlation time. 
This is precisely the decaying time of the correlation functions of the bath. 
If the correlation time of the bath is much shorter than the typical time scale 
for the system dynamics 5 , then we can make in (3.21) the replacement: 

a(t') -> a(t) (3.23) 

and obtain in this way a differential equation for a. Finally, if we are in- 
terested in the dynamics of the system for times much longer than the bath 
correlation time, the lower integration limit in (3.21) can be moved to — oo 
since the initial bath correlation are irrelevant. With this set of approxima- 
tions the knowledge of the state of the system at some time t is enough 
to determine the state at all times t > t . This property is called Markov 
property. In the weak coupling limit and assuming short correlation times 
in the heat bath we have derived the following GME for the reduced density 
operator in the interaction picture: 

1 r°° ~ 

~° = -ft\ drTr B {[#i(i), [Hjit - r), a(t) <g> p B ]]} (3.24) 

To proceed further the precise knowledge of the model Hamiltonian is 
required. In section 3.4 we will specialize this derivation of the GME to the 
description of the dissipative environment of the shuttling devices. 

3.3 Derivation "a la Gurvitz" 

The tunneling coupling of the shuttling devices to their electrical baths (the 
leads) is not weak. It sets, on the contrary, the time scale of the electrical 

4 Note that causality is preserved since the state of the system at times t' > t does not 
enter the integral. 

5 By system dynamics we mean in this case the time evolution of the reduced density 
operator in the interaction picture a(t). In this sense it is the weak coupling to keep the 
system dynamics slow. 
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dynamics that in the shuttling regime is comparable with the period of the 
mechanical oscillations in the system. 

In the SDQS the tunneling amplitudes depend exponentially on the dis- 
placement of the central dot from the equilibrium position. The oscillations of 
the QD modify correspondingly the tunneling rates. In the shuttling regime 
the following non-adiabatic condition is fulfilled: 

r L , R = ^ (3.25) 

where the average can be interpreted as a classical average over the stable 
limit cycle trajectory or quantum mechanically as an expectation value in 
the stationary state. In both cases Coulomb blockade must be taken into 
account for a correct evaluation of the average rate 6 . 

In the TDQS the coupling to the leads is constant and represents a tun- 
able parameter of the model. Also for this device the cleanest shuttling 
regime is achieved for a rather high coupling (and associated tunneling rates 
comparable with the mechanical frequency). 

In 1996 S. Gurvitz and Ya. S. Prager proposed a microscopic derivation 
of the GME 7 for quantum transport with an arbitrary coupling to the leads 
[28]. Following their article we give now in detail the derivation of the rate 
equations for transport of non-interacting spin-less particles through a static 
single dot connected to leads 8 . Even if in this oversimplified case the result is 
intuitive and could be guessed just using common sense, the generic features 
of the derivation will appear. We extend the result to particles with spin 
in strong Coulomb blockade and eventually we conclude the section showing 
that also coherent transport can be treated in this formalism and we derive 
the GME for a static double dot device. 

Let us consider a quantum dot connected to two electronic leads. We 
neglect (at the moment) the spin degree of freedom and Coulomb interaction 
between the electrons. The energy levels in the macroscopic leads are very 
dense while they are discrete in the microscopic QD. We assume that only 
one level in the dot participates in the dynamics of the model 9 . 

6 We will discuss the details in section 6. 

7 In the article they use the expression "rate equations". Nevertheless the equations 
derived fully involve coherencies. 

8 It must be noticed that the problem of calculation of current through arrays of static 
quantum dots has been analyzed in more or less equivalent approach also by other authors. 
We cite as example for similarity of results Wegewijs and Nazarov [29]. 

9 Whilc the non interacting approximation for the leads can be understood in the frame 
of Landau theory of quasi-particles, no interaction combined with single level approxima- 
tion for the QD is an oversimplification with no physical explanation that must (and will) 
be relaxed. 
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3.3.1 Many-body basis expansion 



The Hamiltonian for the model reads 



c, + eicjcj + ^ 



r 



(3.26) 



+ ^fi,(cf Cl + cjc,) + ^fi r ( c t Cl + c J Cr ) 



r 



where cj,q,cj create a particle in the left lead (energy level £/), in the dot 
and in the right lead (energy level e r ) respectively. e\ is the energy level 
of the dot and f^( r ) the tunneling amplitudes between the dot and the left 
(right) lead. For temperatures much smaller than the Fermi energy of the 
leads we can approximate their Fermi distributions by step functions. The 
chemical potential of the left (right) lead is assumed much higher (lower) 
than the dot energy level. 

We identify the empty state |0) for the model with the condition of empty 
dot and leads filled up to their Fermi energies. Then we gradually move 
electrons from the emitter to the dot and finally to the collector and associate 
a new vector to each state of this "decaying chain" . We construct in this way 
an infinite many-body basis that defines the Hilbert space for the model 10 . 
The first few elements of the basis read: 



|0> c\c h \0) 4 lC jO) clcl lCli c h \0) 4A%%\U) ■■■ (3-27) 



where we choose by convention to move all the creation operators to the left 
and the annihilation operators to the right. We also assume with the two 
groups that < Ei i+1 and e Ti < e Ti+1 to avoid double counting. The state 
of the model is described by the many-particle vector which we can 

expand over the basis: 



10 Somc of the possible states of the device are excluded from this Hilbert space: states 
with electrons excited above the Fermi level of the left lead and/or holes below the Fermi 
energy of the right lead. They are neglected because they would anyway be hardly pop- 
ulated due to the fast relaxation of the leads to their thermal states and the very low 
probability of electron tunneling for energies so far from the resonant level of the dot. 



|tt(t)) =[&„(t) + J2 b ^(t)c\c h +$> 1 r 1 (t)cJ 1 C Jl 





(3.28) 
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Figure 3.1: Schematic representation of the first elements of the many-body basis 
for the single dot model. Electrons are progressively taken from the emitter and 
moved to the QD and finally to the collector. The vectors represented are in the 
order a) = |0), b) = c{c h \0), c) = cj 1 q 1 |0) ; d) = cjcJ^qJO), e) = c^c^q.qJO). 
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3.3.2 Recursive equation of motion for the coefficients 

The vector obeys the Schrodinger equation ih\^f(t)) = H\^f(t)) and 

we impose the initial condition |^(0)) = |0). In terms of the coefficients of 
the expansion (3.28) we obtain an infinite set of differential equations with 
the initial condition b (0) = 1 and all the other coefficients equal to at time 
t = 0. 

Due to the quadratic form of the Hamiltonian, the infinite set of differ- 
ential equations for the coefficients 6's presents a recursive structure: each 
coefficient is linked in its equation of motion only to the previous and the 
next in the "decaying chain" (3.28). Since we are interested in keeping track 
of the state of the dot we condense the full system of equations into two 
equations for the generic coefficients bn \n_ s r ,\n_ (t) and b-, ,, .-,„ (t): 



ib {h}u^U 



k=l 



b {lj}^ =1 {rj}^ =1 



in+1 



k = l 



ib 1 



£ l + ^2( £ r k ~ £ h) ~~ £ ln+l 



k=l 



b MUM 



n+1 



(3.29) 



r n +i 



k=l 



where h = 1 and the sums over labels (e.g. X^ n+1 ) are continuous sums 
over all the possible energy levels of the leads. We also used the shortened 
notations: 



b Hh}£t\{lk}{rj}? =1 

In order to proceed in the derivation of the rate equations it is most con- 
venient to make the Laplace transform of the system of differential equations 
(3.29). We obtain the system of algebraic equations: 



I n r 1 r2r :i ...r n { t ) 

(3.30) 

bll 1 l2h...l k - 1 l k+ i...l n +irir2r3...r n \t) 



35 



CHAPTER 3. GENERALIZED MASTER EQUATION 



E + ^2( £ i k ~ £ r k ] 



k=i 



ln+1 



nO 



k=l 



E-e 1 + ^2(e lk -e rk )+e ln+1 



k=i 



n+1 



(3.31; 



r„+i 



k=l 



where the Laplace transformed coefficients are indicated with a tilde and 
are functions of the variable E (that we assume to have an imaginary part 
to ensure convergence of the Laplace integral). At this level the left-right 
asymmetry reveals itself in the number of "decay channels". Each state of 
the chain (3.27) is coupled to an infinite number of right states and only a 
finite number of left states. Since the couplings are equivalent this results in 
a statistically definite direction of motion for the electrons. 



3.3.3 Injection and ejection rates 

The continuous sums in (3.31) can be simplified using the recursive structure 
of the equation of motion. We isolate in the second equation of (3.31) the 
coefficient ^>i{i 3 } n+ l{r 3 } n ( E ) an d insert the result into the first equation of 
(3.31). The continuous sum results into two terms: 



*n+l 'n+1 



and 

n+1 



£7/ I -i f2/ 7 

tn + l 

~ L ' ft] — Ft A- V^™ 

ln+1 k=l 

Since the energy levels in the leads are dense we can substitute 

»+oo 



/+oo 
de ln+1 D L (e ln+1 ) (3.34) 
-oo 

'n+1 
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where D(e ln+1 ) is the density of states in the left lead calculated at energy 
si n+1 and we have extended the integration limits to infinity in the wide band 
and high bias approximation. We can evaluate now the sum over / n +i in 
(3.32) and (3.33) using residues method. Since all the poles are in the same 
half plane we can neglect all terms which are asymptotically o(|e; n+1 | _1 ) for 
\ £ in+i \ ~* 00 • ft i s cl ear from the algebraic system (3.31) that a coefficient b 
that has among its indices £i n+1 behaves asymptotically at least like \ei n 
For this reason (3.32) vanishes and only one term is left from (3.33): 



i-i 



+i i 



2 = 1 



+ 00 



de l+1 D L (e ln+1 ) 



ln+1 



E - ei + Er=i( £ ii - £ n) + £ i n 



(3.35) 



+i 



For the evaluation of the integral is enough that the tunneling amplitude 
Vti n+1 and the density of states DL{ei n+1 ) are analytical and non-zero where 
the denominator vanishes. We assume that they are constant to avoid n- 
dependence in the tunneling rates and perform the integral in (3.35). The 
second equation (3.31) can be treated analogously and the system reads: 



E + ^2( £ h - £ r k )+i^r 



k=i 



J^(-l) fc n n rk b 1{h} n =i{rj} „ =iX{rk} (E) =i5 n0 



k=i 



k=i 

n+l 



(3.36) 



{limi{rM =1 



-E(- 1 ) fc " B " ln '^ft}^Hr J > y=1 (^ = 



k=l 



where we have introduced the injection and ejection rates T L and T R 



T L = 2nD L Q 2 L 
T R = 2tiD r Q 2 r 



(3.37) 



whith the energy independent tunneling amplitudes (Ql and fl R ) and density 
of states (D L and D R ) U . 

11 We assume real tunneling amplitudes as it is also implied by the form of the Hamil- 
tonian (3.26). The most general case of complex amplitude would result anyway in the 
tunneling rates T LtR = 2irD L . R \tl L . R \ 2 . 
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3.3.4 The reduced density matrix 

The reduced density operator is defined as the trace over the bath variables 
of the total density operator: 

a(t) = Tr B {\*(t)}(*(t)\} (3.38) 
The matrix elements of the reduced density operators are explicitly 

= ]Tfe,£W)>W)|fi,Js> (3.39) 

{B} 

where \is, B), % — 0, 1 is the vector that corresponds to the empty or charged 
dot (the system) and a particular configuration B of the leads (the baths). 
We assume that the bath state B does not contain coherent superpositions 
of states with different number of particles. This implies the vanishing of co- 
herencies in the reduced density matrix. It is useful to organize the sum over 
the bath configurations according to the number of extra electrons (holes) 
collected into the right (left) lead. 



MO = EEfeBn|*W>W)|4i S > = 5>£°(*) (3-40) 

"=0 {B n } n=0 

where B n is a configuration of the baths with n extra electrons in the col- 
lector and we have introduced the n-resolved density matrix a^ n \ Using the 
expansion of the vector in the many-body basis (3.27) we can express 

the n-resolved density matrix in terms of the coefficients b. For the two 
non-vanishing elements: 



{B n } {lk}{r k } 

a^ = j2\(±s,B n mt)}\ 2 = E iWjkkJ^ 

{B n } {lk}{r k } 



(3.41) 



where the sums are calculated over all the possible configurations of indis- 
tinguishable particles (e.g. J2{i k } = Zw 1 <i 2 <j 3 < <«„)■ ^ ne ti m e-dependent 
matrix elements of the reduced density matrix are connected to the coeffi- 
cients b(E) by the inverse of the Laplace transform: 



38 



3.3. DERIVATION "A LA GURVITZ" 



{h}{rk} 



>\J(E-E')t 



{l k }{r k }' 



(3.42) 



Apart from being a natural step in the derivation of the GME in the Gurvitz 
approach, the n-resolved density matrix contains the additional information 
on the number of electrons collected in the resevoir at time t. This informa- 
tion is very useful to the calculation of the current noise in the SDQS where 
the quantum regression theorem can not be applied due to the form of the 
current operator that involves both system and bath operators. 



3.3.5 Generalized Master Equation 

The equation of motion for the reduced density matrix is obtained combining 
(3.42) and (3.36). First we derive an equation of motion for the n-resolved 
reduced density matrix a^ n \ The case of the empty dot population with 
n = is special due to the particular choice of the initial condition and we 
treat it separately. The starting point is the first of the equations (3.36) 
specialized for n = 0, namely: 

(E + i^jbo(E)=i (3.43) 
We taking the inverse Laplace transform and obtain 

bo(t) = ~Y b o(t) 

The definition of a^f and the Leibnitz theorem for derivatives 
conclusion: 

4? = b b* + b b* = -T L b b* = -r L a$ (3.45) 

This argument is applied also in the case with n ^ but the struc- 
ture of the equation is more complex and in general a final continuous sum 
must be evaluated. We take the first equation in (3.36) and multiply it 
by by^yn {,..}« ^(E')e~ l(yE ~ E )*. Then we subtract side by side the complex 
conjugate of the first equation of (3.36) evaluated in E 1 and multiplied by 
b{i j } n _ 1 {r j } n _ 1 (E)e~ tlyE ~ E )*. Finally we integrate in dE and dE' and sum over 



(3.44) 
lead to the 
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the bath configurations with n electrons in the collector. We repeat the 
procedure also for the second equation in (3.36) and obtain: 



E [^[0-0 + ir l )6 {ij)?=iH)?=i (£)^ )?=i{rj)?=i (B0e-(— 

n 

-23(£(-l) fe -^ rfc ^^^ 



: 

(3.46) 

for the first equation and similarly 

n+l 

(3.47) 



k=i 



for the second. S indicates the imaginary part. In the definition of the n- 
resolved reduced density matrix the two coefficients b correspond to the same 
bath configuration. The finite sums in equations (3.46) and (3.47) still have 
coefficients with different bath configuration. Using properties of the Laplace 
transform, the definition of and the relations 



ES(-i) fc n_ln ^; i3 . }&luw{rj} , =i (^ 



b* 



(3.48) 



^^ }7 J E ') E ,_ £i + ELi(£ifc _ ^ + _ ^ 
obtained from (3.36), we transform (3.46) and (3.47) into: 



^ + rvC) = 



E 

{lj}{rj} 



dEdE' 
4vr 2 



M y — 

\ £-^< 771/ 



-i) k + k 'n n 



*ift)I =1 h)J =1 \hl( fi )^ ) )» =1 ( r) )? =1 \h4( E ') e " !(,i " ,! ' )t 
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and 

• / • ( n ) i t~i (n) \ 



v [dEd&r , ^ (-l) fc+fc ^A 

(3.50) 

It is crucial at this point that k = k'. If k ^ k! we can eliminate the 
variables (Ik) from 6 and tv (ly) from 6*, perform the integral over one 
of the now common "missing" variables 12 and obtain zero. We are left with 
the case k — k' . We transform the sum over the "missing" variable r\ (Ik) 
into an integral in the corresponding energy. The discrete sum in the index 
k takes care of the integration limits and sets them to infinity. The integral 
can be performed using residues methods to get: 



{IjHrj} 



l [a^(t) + T R a{ 1 l ) (t)]= l T L £ \b m=l{rj}7=1 (t)\ 

Oi }{»-,■} 



2 



(3.51) 



Finally we use the representation of the n-resolved reduced density matrix 
(3.41) and obtain the master equation: 

*» - ~ riff :, + ^l 1 , (3.52) 

where we assume that cr^ ^ = to include into the same compact form also 
the equation (3.45) for n — 0. From this set of equations it is possible to 
determine the current in the left and right leads. The current in the right 
lead is the time derivative of the total number of electrons collected in the 
right lead at time t: 

oo 

I R (t) = N R (t) = $>[4o } W + *&\t)] (3.53) 

?1=0 

12 Missing in the sense that they have been eliminated from the coefficients subcripts. 
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Inserting (3.52) we obtain the intuitive result: 

oo oo 

I R (t) = j^nTfl^r 1 ^*) " <r{?(f)] = r fl $><?(*) = T R a u (t) (3.54) 

n=0 n=0 

For the calculation of the left lead current we have to start with the analog 
of (3.52) but this time resolved for the number of holes h accumulated in the 
emitter: 

h {h) _ r Ah) r (h) 

Jh) _ r Jh) , r _(fc-l) [ } 

a ll — ~ L R a ll + 1 £°00 

The left lead current reads: 

I L (t)=N L (t) = T L a 00 (t). (3.56) 

The average over the bath degrees of freedom is completed by summing 
(3.52) or (3.55) over all the possible number of electrons (holes) collected in 
the right (left) lead. 

o"oo = -TlOoo + r#crn 

(3.57J 

o"n = — r^on + r L a o 



The system (3.57) is a set of rate equations for a two-state model. The 
empty and charged states are connected by charging and discharging rates 
(T^ and respectively) and the variations in the populations of the two 
states is given by a balance of incoming and outgoing currents. The sta- 
tionary solution of (3.57) is achieved for I L = I R 13 . This condition and the 
general sum rule <t o + Cn = 1 give the stationary populations: 



°oo = 



r, 



a 



(3.58) 



11 



r L + r 



R 



and the stationary current: 



I st = 7^^- (3-59) 

13 It is easy to verify that the condition of balanced current is equivalent to the stationary 
condition &qq = &u =0 
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Figure 3.2: Schematic representation of the dynamics of the single dot device 
with spin. The two spin species can tunnel to and from the two degenerate spin 
levels of the quantum dot with rates , T ^ and , , respectively and without 
influencing each other. Due to Coulomb blockade only one electron at a time can 
occupy the dot. 

3.3.6 Spin and strong Coulomb blockade 

The rate equations (3.57) are an intuitive result that can be written sim- 
ply using common sense. Nevertheless the effort spent for their microscopic 
derivation is justified by the possible generalizations that will lead us to the 
GME for shuttle devices. First we want to relax the assumption of spin-less 
non-interacting particles. The spin of the electrons can be very easily taken 
into account if we assume strong Coulomb repulsion in the dot. Due to a 
charging energy much larger than any other energy in the model we assume 
that only one electron at a time can occupy the dot. The Hamiltonian reads: 

H = Y ElC \s C ls + ^As C ls + Y £ r C ts C rs 
l,s r,s 

+ Y ^( c l°is + C UJ + Y n r(ct s c ls + iO (3.60) 

l,s r,s 

+ Uc\ s c ls c\_ s c 1 _ s 

which is an extension of the Hamiltonian (3.26) where s = ±1/2 is the 
spin degree of freedom and U the charging energy of the double occupied 
dot. We take into account the interaction in the definition of the Hilbert 
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Figure 3.3: If we neglect the information of the spin specie the system dynamics 
can be reduced to the one of an effective spinless system with asymmetric injection 
and ejection rates. 

space by discarding all the many-body states with double occupied dot. The 
effective Hamiltonian that we consider is then quadratic, and the Schrodinger 
equation projected onto the many-body basis gives rise to a recursive set of 
equations similar to (3.29). We have in this case three general equations 
corresponding to the three different states of the quantum dot: 



ib 



\Vj} n j l l {lh} n j Utfr j } n j l 1 {lr j } 





(3.61) 




k=l 
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for the empty dot coefficient, 



ib 



k=l 



, „ n-t +1 , , n i nt f - 1 



r n T +l 

n T +l 



k=l 



\{Ti*}{iij}"ii{Trj}"l 1 {lr j }"i 1 



(3.62) 



for the spin-up and finally 



26 



fc=l 



, n t r ni +1 , , nt , 



r n T +l 



+ V(-1)'-Tiaft + 



k=i 



1 \{i'ik}{Tr J -}^I 1 {lr j };i 1 



(3.63) 

for the spin-down coefficient. In the last three differential equations (3.61), 
(3.62) and (3.63) we have extended the notation used in equation (3.29) to 
take into account the spin degree of freedom. Despite the heavy but complete 
notation that keeps track of the four baths (two leads with two spin species 
per lead) and the state of the dot, the same kind of arguments that we used 
for the spin-less case bring us to the set of rate equations: 

°oo - ~\ L + 1 liJVqo + 1 flT°TT + i fli cr u 



(n) 

u 

(n) / (n)x 



(3.64) 



-r R icr 



(n) 



1 M°oo 



where ( (7 }i ) is the population of spin up (down) in the dot with n = 
n I + rij electrons in the collector and we have introduced the spin-dependent 
injection and ejection rates: 
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Fl,k\ = 2-nD LM Vt LR 
Fl,ri = 27rD LiRl Q 2 LR 



(3.65) 



The sum over the number of electrons in the collector gives: 



O"00 



(3.66) 



The coherencies between different spin species in the QD (e.g. a^) vanish 
in this model because different spin states of the dot correspond to differ- 
ent bath states and the only way to have coherent superpositions in the dot 
would be to maintain the same in the leads which instead are assumed (as 
macroscopic objects) incoherent 14 . If we are not interested in the spin in- 
formation on the dot we can introduce the population for the charged dot 
= aff + 0|?\ Assuming also non-polarized leads (i.e. Dl,ri = ^l,ri) 
and, consequently, tunneling rates independent of different spin species the 
system of rate equations (3.66) becomes: 



where F^ = = and Fr = F^ = F^. Comparing these rate equa- 
tions with the ones derived for the spin-less non-interacting model (3.57) we 
note that the only remaining signature of the spin degree of freedom is in the 
injection rate. In the case of identical leads the injection rate doubles the 
ejection rate. This behaviour can be interpreted in terms of tunneling chan- 
nels: both spin species can tunnel in when the dot is empty, but once the dot 
is charged with an electron of specific spin only that species can tunnel out. 
At this level the spin degree of freedom is just renormalizing the injection 
rate. Since this argument can be repeated for any model in strong Coulomb 
blockade we will restrict the derivation of the GME for shuttling devices to 
spin-less non-interacting particles. 

14 Non-trivial spin coherencies can be achieved for example by introducing a spin dy- 
namics in the dot. In that case different spin states on the dot could correspond to the 
same bath state. 



coo = 



2r L (Too + F R a n 
Fr<j u + 2F L a o 



(3.67) 
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Figure 3.4: Schematic representation of the first many-body basis elements for the 
double dot system. 



3.3.7 Coherencies and double-dot model 

A simple example of a device that exhibits coherent transport is represented 
by an array of two quantum dots located between a source and a drain lead. 
We assume that the device is working in strong Coulomb blockade (i.e. only 
one electron at a time can occupy the device, either in the left or in the right 
dot). Electrons can tunnel in the device from the emitter only to the left dot 
while tunneling off is allowed only from the right dot. This condition can be 
achieved due to the fact that the tunneling coupling to the leads decreases 
exponentially with the distance and can be neglected for the far lead. Also 
the two dots are in tunneling contact. Since the transport must happen via 
tunneling between the discrete levels of the dots we expect coherencies to 
play a role. The Hamiltonian for the model reads: 



H =e L \L)(L\ + Q (\L)(R\ + \R)(L\) + e R \R){R\ 

+ V eic]^ + V e r c\c r 

i r (3-68) 

+ J2 n,(|£Xo|^ + \o)(l\4) + J2 Vr(\R)(0\c r + \0){R\4) 

I r 

We identify in the first line the system Hamiltonian with the single energy 
levels of the left and right dot (el and Er) and the tunneling amplitude Qq. 
The second and third lines describe respectively the electronic baths (the 
leads) and their coupling to the device. The system can be found, as in 
the spin model, in three different states: empty or occupied with an electron 
either in the left or right dot. We associate to each of these states the vectors 
|0), \L) and \R). The Schrodinger equation projected onto the many-body 
basis for the system+baths Hilbert space can be then represented by the 
following three recursive differential equations for the expansion coefficients: 
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.fc=l 

+ S^wAft^+V,}^ + ^(- 1 )^ nfi r fc &i?{i J }y =1 {r J }y =1 \K} 
in+l 

n 



k=l 



k=l 



+ Q 0^ {i . } ™+l {r . } n x + $^(- 1 )^' 1_1 ^ fo ft}"+ 1 \{i fc }{r i }^ 1 



n+1 



fc=l 



26 



^■}"±i{^}?=l 



£ R + y^( £ r^ - £« fc ) - £/„ H 



fc=l 



+ Q 0fc Lfe} "+l {r . } n x + ^ fi r n+1 & { i 3} «+i {r:;} "+i 



(3.69) 



The Laplace transform can be taken and the continuous sums in the first and 
third equation be performed to give the set of algebraic equations: 



k=l 



& fc}? =1 {r,}? =1 (£) 

^(-l) k ~ n Q rk b R{h} n =i{rj} n^ {rk} (E) = i5 n0 



k=l 



E + e L + ^2(e rk - E h ) - e ln+1 



k=i 



~ h L{l 3 }-tl{r 3 }^S E ) 



n+1 



n ^ R{lj}1 +i {rj} n J,E) -$^(-l) fe n 1 ^i k b{i j }"+i\{i k }{r j }™(E) = 



k=i 



E + € R + J2( £ rk - £ h) - £ l n+ i + 
k=l 

-n h L{lj} n+i {rj} (E) = 



n (E) 



(3.70) 



In the double dot model some coherencies of the reduced density matrix 
do not vanish since they correspond to different "internal" states of the dot 
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and can be combined with the same state of the baths. For example: 



*U(t) = Y,(L,BnMt)mm,Bn) 



{h}{rk} 

The next step is an n-resolved GME for the reduced density matrix a. 
The equations for the populations are derived following the procedure we 
explained for the single dot model: 

• ( n ) _ T ( n ) i T ( n— 1) 
°00 - _i L&00 + 1 R a RR 

42 = -iOo (42 - 41) + r^o"o (3.72) 

• (™) T-l ( n ) l~> Z' (™) ( n )\ 

^RR = - T R°RR ~ Mo [ a LR ~ RL ) 

The coherencies play an active role in the transport through the quantum 
dots: the second and third equations in (3.72) show that the left (right) dot 
can be discharged (charged) only via coherent transport. We concentrate 
now on the equation for the coherence 41- We take the second equation in 
(3.70) and multiply it by b* xn+1} xn (E')e-^ E - E '^\ then we subtract side 

by side the complex conjugate of the third equation in (3.70) evaluated in E' 
and multiplied by ^L{i j } n+ l{r j } r !- i (E)e~^ E ~ E '^ t . Finally we integrate over dE 
and dE' and sum over the baths configurations with n electrons in the col- 
lector. Using the properties of the Laplace transform and the representation 
of the reduced density matrix in terms of the coefficients b of the many-body 
expansion we obtain: 

r, 



+ (*l -e R + i-fytkt) - n [41(0 - 42 w 

{l J }{r J } J k=l 



The last term in the LHS of equation (3.73) vanishes since the integrand 
behaves asymptotically as o{\ei., r . | _1 ) in the limit |£/-, r -| — > oo for all the 
integrated variables £j. >r ... The average over the bath degrees of freedom is 
completed by the sum over the number of electrons in the collector n that 
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leads to the GME: 



o"oo = — Tl^oo + F R a RR 
oll = -i^oivRL - cr LR ) + T L a 00 

&RR = -TrVRR - iQo{&LR - CTRL) 

o'er = -i(£L - £r)clr - |r R a LjR - iil (a RR - a LL ) 

where we have omitted the equation for ct R l since a is a hermitian operator 
on the Hilbert space of the system and <jrl = &lr- We can better visu- 
alize the coherent and incoherent contribution to the GME using a matrix 
representation: 

& = -i[H sys , a] + E[a] (3.75) 
where a is the density matrix 

°oo °ol cor \ 

o-lo o- LL a LR , (3.76) 
pro & rl &rr ) 

H sys is the Hamiltonian for the system extended to the empty state for the 
system 

/ \ 

H sys = e L n \ (3.77) 

V o n e R J 

and S[«] is a linear super-operator that transform operators on the Hilbert 
space of the system into operators on the same space and acts on the density 
matrix: 

-r L o- 00 + T R a RR 
3[<t] = ( IVoo -\V R a LR \ (3.78) 

— -^R^RL —F R a RR 



3.4 GME for shuttle devices 

The shuttle devices are in contact with two different kinds of bath: two 
electrical baths (the leads) and a mechanical bath. We assume that the 
electrical and mechanical baths act independently on the device. This as- 
sumption splits the GME into two additive components, one for each kind of 
bath. Due to the different coupling strengths we derive the electrical compo- 
nent extending the method proposed by Gurvitz and extensively presented in 
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Section 3.3 while for the mechanical component we adopt the weak coupling 
quantum optical derivation presented in Section 3.2. 



3.4.1 Single Dot Quantum Shuttle 

We start recalling the Hamiltonian for the SDQS: 



H — H sys + H\ eads + ifbath + H tnn + H\ 



1 int 



(3.79) 



where 



(3.80) 



p 2 1 

-ffsvs = t; 1- -mu 2 x 2 + (ei — e£x)c\c-i 

2m 2 

Pleads = ^{£lA k C l k + £ rA k C r k ) 
k 

H tun = J2[Ti(x)cl kCl + T r (x)ct kCl ] + h.c. 

k 

#bath + H int = generic heat bath 



Also the mechanical degree of freedom is treated quantum mechanically. 
For example the position operator can be expressed in the form: 



x 



h 



{£ + d) 



(3.81) 



where S and d are respectively the creation and annihilation operators for 
the harmonic oscillator. We neglect for a moment the mechanical bath and 
its coupling to the system and start the Gurvitz analysis of the model dy- 
namics. The many-body basis introduced in (3.27) must be extended to take 
into account also the phononic excitations of the system. For definiteness 
we choose the eigenvectors of the oscillator Hamiltonian as a basis for the 
mechanical part. We display the basis elements in the following table: 



r=0 



|0> 

c\c h \0) 

4^10) 



dt|0) 



S\o) 



c i c h 
4^/10) 



r=2 



<no) 



^10) 



(3.82) 



where r is the number of excitations of the oscillator and the empty state 
|0) represents the empty dot in its mechanical ground state with leads filled 
up to their Fermi energies. It is convenient to organize the coefficients of the 
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expansion of the state vector |\&) in the basis (3.82) in vectors: one for each 
electronic configuration. The different elements of the vectors refer to the 
different excited states for the oscillator. The Schrodinger equation for the 
state vector is represented in the basis (3.82) by two recursive differential 
equations for the vector coefficients b's: 



ih 



{'j}?=ifo}S=i 



-^osc + ^ ] ( g r fc 



k=l 



tn+l 



lb, in 



+ ^2 T l(x) h l{l j }]+l{r j }J =1 + ^(-^^ " T r(^)bl {ij .}n =i{r . } n =iXK} 

k=l 
n 

H osc + ei - e£x + ^(£r fe - e lk ) - e ln+ 



k=l 
n+1 



'l{Ii}S{^i}" = l 



+ ^T r (£)b fe} , + i {rj} , + i +^(-l) fc n 1 ^(^)b . }? +i\ {ifc}{r . } n =i 

(3.83) 



r n +i 



k=l 



where x is given in its matrix representation in terms of the occupation 
number basis (h — 1): 



2mcu 



(8 r ,s+l + ^r,s-l) 



(3.84) 



and the Hamiltonian for the harmonic oscillator in the same basis reads 15 : 



(3.85) 



All the constants in equation (3.83) are identity operators in the mechan- 
ical Hilbert space. 

One of the key assumptions in the derivation of the GME in the Gurvitz 
approach is the position of the energy levels of the system: they must lie 
well inside the transport window open between the chemical potentials of 
the leads. Since the oscillator spectrum is not bounded from above we as- 
sume that only a finite number of mechanical excitations are involved in 
the dynamics of the system. We will see that, at least in the presence of 
a mechanical bath, this assumption is numerically fulfilled. In any case a 



15 The representation given in equation (3.83) is actually independent of the basis for 
the oscillator Hilbert space. The b vectors are projections of the state vector |\I>} on the 
particular subspace given by the electronic configuration specified by the subscript. 
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violation of this condition in the final result would be unacceptable since it 
would violate the validity condition of the GME. From the point of view of 
the experimental realization of the device this limit is imposed at least by 
the leads that set an upper bound to the amplitude of the dot oscillations. 

The Laplace transform of (3.83) with the initial condition \ty(t = 0)) = |0) 
reads: 



E + H osc + ^2(e rk - si k ) 



k=i 



(n+1 



k=l 



E + H osc + ei - e£x + ^2(e rk - e tk ) - e tn+1 



k=i 



n+1 



(3.86) 



r n +l 



k=l 



where vo is the initial condition of the oscillator. 

The continuous sums in the system (3.86) can be performed with an 
argument similar to the one used for the static QD. We have just to be 
careful with the matrix notation and change the basis to diagonalize the 
matrix: 



M 



E + H osc + si - e£x + ^2(s rk - Ei k ) - Ei n+1 



k=l 



(3.87) 



before taking the integral. The sum ^ r in the second Eq. of (3.86) can be 
treated analogously. As in the static case the continuous sums are condensed 
into "rates" 16 : 



16 These "rates" are position dependent and then in our quantum treatment they are 
operators. Actual rates can be recovered by averaging these operators on a given the 
quantum state. 
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E + H osc + J2( £r * ~ £l ^> ~ i 



T L (x) 



k=i 



^{hn=i{rjn =1 ( E ) 



£(-l) fe n T r (x)h 1{lj}u{rj}u \ {r . k} (E) =i5 n0 v 



k=i 



E + H osc + ei - e£x + ^2(e rk - 



k=l 



n+1 



k=l 



where 



F LtR {x) = 2irD LtR T? r (x) = T LtR e^ 2£/x 



(3.88) 



(3.89) 



and we have introduced the tunneling length A and the bare injection (ejec- 
tion) rate T L (T R ). 

The reduced density matrix for the system contains information about 
the electrical occupation of the QD and its mechanical state. Coherencies 
between occupied and empty state vanish because they imply coherencies 
between states with different particle number in the baths. The equation 
(3.40) for the non- vanishing elements written in the static QD case still holds: 



(3.90) 



with the difference that the "elements" an are now full operators in the 
mechanical Hilbert space. It is useful to express them in terms of the vectors 
b: 



= £ b {'*}2=i{r*}£ =1 (*)b{i t }» =1 {r t }5J =1 (*) 



{lk}{rk} 



{h}{rk} 



(3.91) 



The notation of equation (3.91) can be understood in terms of the Dirac 
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notation: b+ is the bra of the corresponding vector b (the ket) 17 . The inverse 
Laplace transform brings us back to the vectors b: 

«£?(*> = E b f w tl(rt)J=1 (£)b; W!= , h)tl (£')e-<—')' 

{h}{rk} 

«#<«> = /t^ E K H::l M : J*! wgl *(^ E '" 



(3.92) 

The case with n = must be treated separately. The inverse Laplace 
transform of the first equation in the system (3.88) specialized for n = 
reads: 

ib = H osc b - i^^bo (3.93) 



and its Hermitian conjugate: 

-ibl = blH osc + ibl 1 ^- (3.94) 



where we have used the property of adjoint of vectors and operators (AB)^ = 
B^A^ and the fact that the oscillator Hamiltonian and position operator are 
Hermitian on the mechanical Hilbert space. The combination of equations 
(3.91), (3.93), (3.94) and the Leibnitz rule for derivatives extended to the 
tensor product between vectors and linear forms lead to the first component 
of the GME for the SDQS: 

*$S(t) = -i[H osc ,*i$] - ^f{e~ 2£/ \^} (3.95) 

where [A, B] = AB - BA is the commutator and {A, B} = AB + BA the 
anticommutator between the operators A and B. Equation (3.95) already 
contains the essence of the driving part of the GME: a coherent evolution 
represented by the commutator with the oscillator Hamiltonian and a non- 
coherent term due to the interaction with the bath. In this second con- 
tribution the quantum features are given by the particular ordering of the 
operators. 

For the general case with n^O the procedure is to take the first equation 

of (3.88) evaluated in E and make the tensor product with bt, xn s xn (E') 

\h ij=i vj/j=i 

17 In other terms the linear operator a u is the tensor product of the vector b and the 
linear form . 
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e %( - E E '^, subtract side by side the adjoint of the first equation in (3.88) 
evaluated in E' multiplied (from the right) with the vector ^{i J }^_ 1 {r j }^_ 1 {E) 
e -i(E-E')t^ i n t e g ra te in dE and dE' and sum over the possible bath configu- 
ration with n extra electrons in the right lead. Using the properties of the 
Laplace transform and the representation of the reduced density matrix in 
terms of the vectors b (3.91) we get: 



lu 00 ^030^00 J 2 I / 



dEdE'* 

4tt 2 



^2(~l) k n r r(^)bi { i.}" =l{r . } n =iUrfc} (S)b { ,. } „ ^S') 



k=l 



-E 

{IjHrj}' 

- t ft}^{ri^(^^}J =1 {r J }? 1 \h}( £ ') r '-( £ ) 



e -i(E-E>)t = o 



(3.96) 



We solve the first equation in (3.88) with respect to b{; i }«_ 1 {r J }^_ 1 - Then 
we insert the result and its adjoint in (3.96) and, as in the static QD, we are 
left with the only non-vanishing continuous sum in the "missing" variable r n . 
The result is the matrix equation: 



- ^{e- 2£ /\^} + r fle */V&-V/* (3-97) 

The treatment of the second equation in (3.88) is totally analogous and 
brings us to an equation of motion for the charged component of the density 
matrix (cr^). Collecting all the results we can write the n-resolved GME: 



°00 — 



1J -OSC) U QQ 



°00 



°11 



—I 



—I 



n °00 



- T t{ 



• * ^' ( hi } • 1'/," " ■ 



Hose cSx, 



-2x/\ ^(n) 



+2x/X 



+ T L e 



i/xM £ /x 



(3.98) 



°oo e 



In order to complete the description of the dynamics of the SDQS we 
have to take into account the mechanical bath and its interaction with the 
system. We derive the mechanical component of the GME starting from 
the general formulation for the equation of motion of the reduced density 
matrix (3.24). We consider the problem described by the Hamiltonian (at 
the moment independent of the electronic dynamics): 



H = H sys + H] 



'bath 



hit 



(3.99) 
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H svs = fojj ( \ + d)d 



where 

'sys \ 2 

7 (3.100) 

#bath = ZlCJqrfq^q 

q 

and the generic interaction contribution: 

H int = hJ2G a F a (3.101) 

a 

G a being a system operator and F a a bath operator and a a generic quan- 
tum number. We will specialize later the interaction in the form we have 
introduced in the chapter dedicated to the model: 

#int = + d\){d + S) (3.102) 

q 

and with the rotating wave approximation 

HlT A) =^9«d + d %) (3-103) 
q 

We start recalling the GME (3.24) derived in the weak coupling using the 
quantum optical formalism: 

i r°° 

~° = -ft\ drTr B {[H int (t), [H int {t - r), a(t) ® p B ]]} (3.104) 
With the generic form for the interaction Hamiltonian (3.101) we get: 

POD 

Ht) = ~ / rfrV{[G a (t)G b (t-r)a(t)-G 6 (t-r)^(t)G a (t)](F a (r)F 6 (0)) 

+[a(t)G a (t - r)G b (t) - G b (t)a(t)G a (t - r)](F a (0)F b (r))} 

(3.105) 

where the tilde indicates the interaction picture and (•) = Ttb{pb*}- We 
can easily go to the Schrodinger picture: 

a = --[H sys ,a} + V / dr{[G 6 (-r)(F o (r)F b (0))(T,G a ] 

TT^' (3.106) 

+ [G 6 ,aG a (-r)(F a (r)F 6 (0))]} 
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Following [30] we introduce the compact notation: 



a = ~[H sys ,a] + £{[G+<7,G«] + [G a ,aG-]} (3.107) 



where 



POO 

G+ a=Y, drG b (-T)(F a (r)F b (0)) 
b Jo 

poo 

G a=J2 drG b (-T)(F b (0)F a (r)) 
i. Jo 



(3.108) 



This formalism is very efficient since, for a given interaction, it requires for 
the GME simply the calculation of the two operators G + and G~ . For the 
interaction Hamiltonian (3.102) we identify 



G q = d + 

and calculate 



(3.109) 



_ POO 

G q=E/ ^ q ,(-T)(F q (T)F q ,(0)> 

q' J ° 

poo 

=d g 2 J dre^ie^riB^) + e~^[l + n B (u; q )]} (3 nQ) 

POO 

+ Sg 2 / rfre-^ r {e^ T n B (cj q ) + e -^ T [l + n B (cu q )]} 



=7T5f 2 {cf[l + n B (c^ q )] + d t nB(o; q )}J(a; - cu q ) 



where we assumed the bath in thermal equilibrium and calculated the aver- 
age: 



* q )=^W-^— y (3.1H) 
and we integrated the exponentials 

jf dre i(lv±UJ ^ T = 7i5(u ±u (l )+V ^ ) ~ ™*(^ ± ^q) ( 3 - 112 ) 

where P denotes the principal value. We neglected the contribution due the 
principal value that only slightly shifts the oscillator frequency oo. Terms 
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proportional to 8{uj + c<j q ) vanish since both uj and w q are positive. The 
operator G~ can be calculated in a similar way and reads: 

G q = V^ubK) + <jft[l + n B K)]}5(w - cu q ) = (G+y (3.113) 
Substituting all the G operators, the GME (3.107) takes the form: 

a = —-[H sys , a] + -[1 + n B (uj)}[d + d\ ad'' - da] + -n B (oj)[d + d\ ad - d)a] 

it Zi Zi 

(3.114) 

where 7 = 27rD(u)g 2 is the damping rate and D(u) is the density of states of 
the phonon bath at the frequency of the oscillator. We rewrite the previous 
equation in the form: 

a = C[a] = C coh [a] + £ damp [ff] (3.115) 

The linear (super) operator C is also known as the Liouvillean and is a linear 
operator defined on the space of linear operators on the Hilbert space of the 
system. The first term £ co h = — ^[H sys ,a] describes the coherent dynamics 
of the isolated system. The terms proportional to 7 and grouped in £dampC 
represent the interaction with the bath which is damping the oscillator. This 
interaction introduces decoherence in the system in the sense that no matter 
how we prepare the initial quantum state of the oscillator (a(t = 0)), in 
absence of other driving forces, the stationary state reached at long times 
(a(t = 00)) is a thermal distribution that corresponds to a diagonal density 
matrix with no coherencies left. The thermal stationary solution is a typical 
property required from a generalized master equation. The GME (3.114) is 
also translationally invariant as can be more directly checked introducing the 
position and momentum operators for the oscillator x,p 18 : 



n B (uj) + 



[x,[x,a\\ (3-116) 



Unfortunately a translationally invariant GME with a thermal stationary 
solution is generating density matrices which are not a priori always positive 
definite. This is a general problem of the GME [30]. In our specific case 
though we checked numerically that in the relevant cases the positivity was 
not broken within numerical errors. 



18 From now on we drop for simplicity the hat for the operators. It will be clear from 
the context if we are dealing with operators or with classical variables. 
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The interaction Hamiltonian in the rotating wave approximation (3.103) 
can be treated in a similar way. One has to extend the space of quantum 
numbers and define: 



G m =d 



Gq 2 = d ] 
F ni = gd\ F q2 = gd n 

The corresponding G + and G~ operators read: 

G+ = dng 2 [l + n B (cu q )]5(cj - cu q ) 

We insert the operators in the general GME (3.107) and obtain: 



(3.117) 



(3.118) 



& = - l -[H sys ,a] + ^n B (u)([Sa,d] + [d\ad]) 

+ ^[l + n B (Lu)]([d,aJ] + [da,J}) 



(3.119) 



where we have defined as usual the damping rate 7 = 1-nDg 2 . Collecting all 
the terms we can write the GME for the SDQS in the form: 



o"oo 



0"11 



El 

2 



[e ,000} +r jR eA(T 11 eA + C damp [cr 00 ] 



(J 11 



e a ,(7n 



^a 00 e a + £ damp [<7n] 



(3.120) 

where we have taken the sum over the number of electrons collected in the 
resevoir and we have introduced the generic damping Liouvillean £damp- One 
can use for example one of the two damping Liouvillean we have just derived. 
In the rest of the thesis we will always adopt for the SDQS the translationally 
invariant damping Liouvillean (3.114). We will also refer to the electronic 
part of the equation (3.120) as to the driving term. In a compact form: 



& = C coil [a] + £driv[o-] + £damp[c] (3.121) 

where we have introduced a block matrix structure to take care of the me- 
chanical and electrical degrees of freedom simultaneously. 
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3.4.2 Triple Dot Quantum Shuttle 

The driving term of the Liouvillean operator for the Triple Dot Quantum 
Shuttle can be derived in strict analogy with the fixed double dot. The 
major simplification with respect to the SDQS is the drop of the position 
dependence in the coupling to the leads as one can see from the Hamiltonian 
for the model: 

H = H sys + H\ cads + -f^bath + H tun + H int (3.122) 

where 

H sys =e |0)(0| + ^f\L)(L\ - ^\C)(C\ - ^\R)(R\ + fiw (dU+^j 
+ t R (x)(\C)(R\ + \R)(C\)+t L (x)(\C)(L\ + \L)(C\) 

Pleads = J2^ C l°l k + e -A C r k ) 
k 

q 

//tun =5>(cJj°X£| +c lk \L){0\) +T r (cl k \0)(R\ +c rk \R){0\) 
k 

q 

(3.123) 

The reduced density matrix for the triple dot system takes into account 
mechanical and electrical degrees of freedom. As in the case of the fixed 
double dot we can organize the density matrix in a block representation: 

o"oo °ol o"oc °OR 

&L0 &LL &LC ®LR 

o = 

CCO OCL &CC &CR 

&R0 crl crc crr 

Due to the incoherent leads the elements cr 0a and a a0 with a = L,C,R 
identically vanish 19 . In the same matrix representation we write 20 the driving 
Liouvillean: 

19 The incoherent leads do not have coherent superposition of states with different parti- 
cle number. Only this kind of "forbidden" bath states would allow coherent superposition 
of states with and 1 particle in the device and thus non-vanishing coherencies a a o or 

20 This equation was first derived in the Gurvitz scheme by Armour and MacKinnon in 
[7]- 



(3.124) 
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--driv 



FrCTrr — T L a 00 








TlOoo 



h^RCTRL —t;Tr(Jrc 





-\^R&LR 
-\^R&CR 
-Fr<J RR 



(3.125) 



where the injection and ejection rates have the form: 



r L>R = 2nD LtR T? r 



(3.126) 



The overall structure of the driving component of the Liouvillean can 
be understood in terms of "decaying channels", since the interaction with 
the continuum of states in the leads gives rise to incoherent tunneling pro- 
cesses. This concept is underlying the following formulation of the incoherent 
dynamics [28]: 



(3.127) 



a'(3'^af3 



where a, a', (3, /3', 5 = 0,L,C,R and T a ^p is the probability per unit time 
for the system to make a transition from state \ct) to state \f3). The generic 
state \a){(3\ is emptied (first line in equation 3.127) and pumped (second line) 
with different T rates. The central dot does not contribute to this incoherent 
dynamics since it is coupled only to the left and right dot discrete states. Due 
to the left-right asymmetry only the rates T ^l = Tl and T R ^ = T R are 
non-zero. The mechanical state of the system is not involved in the equation 
(3.127) but plays an active role in the coherent dynamics. In block matrix 
notation the system Hamiltonian takes the form: 



H S y S 







4p + H osc 







t L (x) 






t R (x) 

AV i TT 

2 ~r ^osc 



t R (x) 
The corresponding coherent Liouvillean reads: 

( £coh[c] ) = — i 2, [(Hsys)aSO'8l3 ~ 0- a s(H sys ) S/3 ] 
V /a/3 „ 



(3.128) 



(3.129) 
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We assume for the damping term the same used by Armour and MacKinnon. 
It is the standard quantum optical damping in RWA that we derived in the 
previous section (see eq. (3.119)): 

£dam P [°"] = — ^rnB(oL>)(dd)a — 2d) ad + add)) 

(3.130) 

- ^[n B {uj) + l]{dUa - 2doS + oSd) 

We write for completeness the system of equations for the 10 non-vanishing 
sub-matrices of the reduced density matrix: 



^00 — —i[H osc , O"0o] + FrCT RR — T L a 00 + £damp[coo] 

o'll = -i[H osc , cr LL ] - it L (x)a CL + ia LC t L (x) + T L a 00 + £ damp [a LL ] 



&cc = —i 



H osc - °cc\ ~ it L (x)cr LC + icr CL t L (x) 



-it R (x)a RC + ia CR t R (x) + C 

damp 

vrr] — it R (x)<7 CR + icr R ct R (x) — T R a RR + £dam P [&rr] 
<?lc = -i[H osc , a LC ] - i^faLC (l + - it L (x)a C c + ^liIl^) 

+iVLRtR(x) + £damp[0"Lc] 

<?cl = -i[H osc , ct C l\ + i^f ( l + ^) °cl + i(?cctL{x) - it L (x)a LL 
-it R (x)a RL + C damp [a C L] 

<?lr = -i[H osc , a LR ] - iAVa LR - it L (x)a C R + i°LctR{x) 
-^a L R + C damp [a LR ] 

<?rl = -i{H osc , <r RL ] + iAVa RL + icr RC t L (x) - it R (x)a CL 

°cr = -i[H OS c, ct C r] - it L (x)a LR - i^f (l - o C r 
-it R (x)a RR + ia C ctR(x) + £ d am P [&cr] 

<?rc = -i[H osc , °rc] + ioRLthix) + i^ORC (l - ^) 

+iv RR t R (x) - it R (x)a C C + £damp [<7rc\ 

(3T3 
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3.5 The stationary solution: a numerical chal- 
lenge 

The master equation generally describes the irreversible dynamics due to the 
coupling between the system and the infinite number of degrees of freedom 
of the environment. It is reasonable to require that in absence of a driving 
mechanism the system tends asymptotically to thermalize with the bath. 
This property is reflected in the evolution of the reduced density matrix 
that in this case has a thermal stationary solution. In the case of shuttle 
devices the oscillator is driven by the electrical dynamics: every time an 
electron jumps onto the moving island it feels an electrostatic force that 
excites the oscillator. For this reason, at least for small enough damping we 
do not expect the oscillator to relax to the stationary thermal distribution. 
Nevertheless since both the electronic and the mechanical degrees of freedom 
of the system are coupled to baths we do expect a stationary solution for the 
GMEs (3.120) and (3.131), i.e. a matrix a stat that fulfills the condition: 



3.5.1 A matter of matrix sizes 

We calculate the stationary matrix numerically: we have to find the null 
vector of the matrix representation for the Liouvillean super-operator C 
The challenge arises from the matrix sizes. In principle the reduced density 
matrix has infinite size due to the mechanical degree of freedom. In order to 
treat the problem numerically we truncate the corresponding Hilbert space 
retaining only the first N states of the harmonic oscillator basis 21 . The size of 
the reduced density matrix a is in the SDQS and TDQS respectively 2N x 2N 
and 4N x 4A: prefactors 2 and 4 come from the size of the electrical Hilbert 
space which is spanned by the vectors |0), |1) for the single-dot device and 
|0), \L), |C), \R) for the triple dot. 

The Liouvillean is a linear operator on the space of Hilbert-Schmidt oper- 
ators 22 on the Hilbert space for the system (Liouville space). Equipped with 
the scalar product: 



21 The case of a different mechanical potential is not more difficult in principle, once we 
know the eigenvectors and eigenvalues of the corresponding one-dimensional Schrodinger 
equation. 



22 Given some Hilbert space H an operator A is of Hilbert-Schmidt if it is linear and 
Tr{AtA} is finite. 



Ca stat = q 



(3.132) 



(A,B) = Tr{A^B} 
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the Liouville space becomes a Hilbert space. One can therefore introduce an 
orthonormal basis and represent linear operators as matrices. The truncated 
Hilbert space for the system gives rise to finite size Liouvillean matrices: for 
the SDQS we reach the size of 4N 2 x 4N 2 while the richer electronic structure 
of the TDQS is reflected in a 16N 2 x 16N 2 elements Liouvillean. Even with 
the condition of incoherent baths that prevents coherencies within states with 
different electron number in the system and therefore sets to zero all sub- 
matrices in the form <r i or a w in the SDQS and a 0a or a a0 with a = L,C, R 
in the TDQS we can not reduce the size of the Liouvillean matrix to more 
than 2N 2 x 2N 2 and KW 2 x KW 2 respectively. 

The description of the shuttle device dynamics requires (especially in the 
shuttling regime) amplitude oscillations of the vibrating dot between 5 and 10 
times larger than the zero point fluctuations. For this reason, in both devices, 
we are left to study the null space of matrices of typical size of 2 • 10 4 x 2 • 10 4 . 
To indicate that we are treating the truncated matrix representation of the 
original stationary state problem we change slightly notation and formulate 
the numerical problem: 



Lp stat = (3.134) 

with Tr{p s * a *} = 1. The solution to this numerical problem came from 
prof. Timo Eirola, Helsinki University of Technology, in the form of a sug- 
gestion and implementation for the SDQS of the iterative Arnoldi scheme. 
We successfully extended the method to the TDQS. 

The Arnoldi scheme applied to the calculation of the null space has ad- 
vantages with respect to the singular value decomposition both in terms of 
computational speed and memory consumption 23 . First we do not need to 
store the Liouvillean matrix and we can always work with operators on the 
system Hilbert space only; second we iteratively look for the best approxima- 
tion to the null vector p stat in spaces which are typically much smaller that 
the Liouville space. Good introductions to the Arnoldi scheme can be found 
in different places in the literature (for example [32] or [31]). We dedicate the 
next section to a detailed analysis of the Arnoldi iteration scheme referring 
for examples to the SDQS Liouvillean. Some of the material can be found 
also in the appendix of [33]. 



23 For the theory and implementation of the Arnoldi scheme we refer to the lecture notes 
"Numerical Linear Algebra; Iterative Methods" by Eirola and Nevanlinna [31]. 
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3.5.2 The Arnoldi scheme 

The central role in the Arnoldi scheme is played by Krylov spaces. For a 
given Liouvillean L and a vector x of the Liouville space we define the 
Krylov space as: 

/Cj(L, x ) = span(x , Lx , . . . , L- 7 ~ 1 x ) (3.135) 

where j is a small 24 natural number. It is important to note that for the 
construction of the Krylov space all what we need are the vectors xo, Lxo, 
L 2 xo, . . . and not explicitly the matrix L. In the SDQS device we would take 
an arbitrary state represented by the two matrices o"oo and o\\ and choose 
a reshaping procedure to map them into a sing le vector x e C 2N ' xl . The 
vector Lx is obtained applying the operator defined in equation (3.120) to 
the density matrices a 0Q and a"n and then reshaping with the same procedure 
used for x . 

The Arnoldi iteration starts with the construction of an orthonormal basis 
{ c lA:}i=i f° r the Krylov space using standard Gram-Schmidt orthonormaliza- 
tion: 



x 
ll x o|| 

q fe+1 = Lqfc ~ E r [q |" Lqfc]qt . k = 1, . . . , j - 1 (3 ' 136) 

ll L qfe-Ei=i[qJ -Lqfclqill 

where the norm in the Liouville space is defined from the canonical Hermitian 
product ||a|| 2 = V at a and has nothing to do with the scalar product we 
introduced to demonstrate that the Liouville space is a Hilbert space. Each 
new tentative basis vector Lq fc is first orthogonalized (by subtracting all 
components in the preceding vectors of the basis) and then normalized. This 
requires the calculation of the quantities: 

hi,k = q\-Lq k ,i = 1, . . . , k; k = 1, . . . , j (3.137) 

and 

k 

h k +i,k = ||Lq fe - 5^[ql • Lq fc ]q;|l> k = (3.138) 

i=l 

that are stored into the upper Hessenberg matrix 



24 We mean small compared to the dimension of the Liouville space. We used j = 20 for 
a Liouville space dimension of roughly 2 • 10 4 . 
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H, 



h 



1,1 


hi,2 


h,3 ■ ■ 


■ Ki 


'2,1 


h 2 ,2 


h 2 ,3 ■ ■ 


■ h 2,j 





h,2 


h,3 ■ ■ 










h 4 ,3 ■ ■ 


h 4 j 



G C j+lxj 



(3.139) 



••• h j+hj 
while the basis vectors are stored as columns in the matrix 

Q, = [ qi ,q 2 ,...,q,]GC 2 ^ (3.140) 

which fulfills Q^Qj = I, I G C jxj being the identity matrix, since the basis 
is orthonormal. The method proceeds by looking for the best approximation 
of the null vector for the Liouvillian within the Krylov space /Cj(L,x ). We 
call this vector Xj. In terms of its j components in the orthonormal basis 
Xj = QjVj. The j coordinates in the Krylov space Vj solve the minimum 
problem: 



min ||Lx|| 2 = ||Lxj|| 2 = ||LQ-Vj||2 (3.141) 

x G Kj 

\\X\\ 2 = 1 

In the process of finding these coordinates a key role is played by the Hes- 
senberg matrix since the following relation holds: 

LQ^- = Qj+iHj (3.142) 

We refer to the notes by Eirola [31] for a rigorous mathematical proof of 
the relation (3.142) and we only limit ourselves to exploit here some of its 
consequences: 



min ||Lx|| 2 = HLQ^Vj^ = ||Qj+iHjVj|| 2 

x G Kj 

l|x||2 = 1 



= v (Qi+iHiVj) 1 • (Qi+i H i v i) 

(3.143) 



(H, Vj )t(Q] +1 Q J+1 )(H,v,) 
J (HjVjO^CHjVj) = \\HjVjh = min ||H jV || 2 



V6P 
IM|2=1 



We started with a minimum problem involving the vector x of length 2A^ 2 
and a matrix L of size 2N 2 x 2A^ 2 and we have reduced it to the minimum 
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problem in the last line that only involves a vector v of length j and a 
matrix H,, of size (j + 1) x j. Since j = 20 the latter is absolutely not 
demanding neither from a computational or a memory point of view on a 
modern computer. We proceed to the minimization using singular value 
decomposition (SVD) [31]. Given a complex matrix A G C mxn with m > 
n, there exist two unitary matrices U G C mxm and V nxn such that A = 
UEV 1 ", and £ G M mxn is diagonal with non- negative diagonal elements 
(conventionally in decreasing order starting with the highest singular value 
in the upper corner [32]). Thus for the Hessenberg matrix H,, 



H, = USV f 



(3.144) 



where U G C^ x(j+1 \ V G C jxj and 



£ = 



Ai 








A, 




G R^ 1 )^' 



(3.145) 



contains the singular values Xj > 0. The minimum problem is solved as 
follows: 



min ||H,-v|| 2 = min HUEVNr^ = min HEVV^ = \ (3.146) 

vec v ea ' v e 

||v|| 2 =l l|v||2=l l|v|| 2 =l 

where we have used the unitarity of U and V and we have chosen the mini- 
mizing vector Vj to be the column of V corresponding to the smallest singular 
value. Having found Vj the best approximation of the null vector within the 
Krylov space can be calculated Xj = QjVj. Finally we test the result and 
compare ||Lxj|| 2 with a given tolerance. If the test is positive we accept the 
result and reshape the vector x 3 - into the matrix form as an approximation 
of the stationary solution of the GME. Otherwise we use x^ as the starting 
vector for a new iteration of the Arnoldi scheme. We have chosen as tolerance 
parameter l(k||L|| 2 where e is the machine precision and the norm 25 of the 
Liouvillean was estimated by T. Eirola as ||L|| 2 ~ exp(A r /log N). 

3.5.3 Preconditioning 

The Arnoldi scheme is iterative and can suffer from convergence problems. 
It is not a priori clear how many iterations one needs to converge and fulfill 

25 We used for the Liouvillean the norm ||L|| 2 = max|| x || 2=1 ^nnrp 
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the convergence criterion. A possible answer to a non-convergent code is to 
reformulate the problem into an equivalent and (hopefully) convergent form. 
This was exactly the situation we had to face with the Arnoldi scheme applied 
to the problem of finding the stationary reduced density matrix for shuttle 
devices. The solution to this problem came again from prof. T. Eirola in the 
form of a preconditioner. The basic idea is to find a regular operator 26 M 
on the Liouville space, invertible, easy to implement, such that the original 
problem C[a stat ] = can be recast into the form: 

M[C[a stat \\ = (3.147) 

and that the finite version of the operator M.C gives rise to a (fast) convergent 
iteration scheme. The operator Ai is also known as the preconditioner. 

The Arnoldi scheme is particularly efficient in finding the best approx- 
imation of the eigenvalues and corresponding eigenvectors for those eigen- 
values that are separated from the rest of the spectrum. Since we want to 
calculate the null vector it is important that the preconditioner moves the 
non- vanishing part of the spectrum far from the origin 27 . 

In order to understand the preconditioner that we used for the shuttling 
problem, we introduce the generic Sylvester operator : C nxm — > C nxm : 



0(X) = AX - XB (3.148) 

where A e C nxn , B e C mxm , X G C nxm . This operator is invertible if and 
only if the spectrua of A and B have empty intersection and the inversion 
routine is computationally light. A part of the Liouvillean operator is of 
Sylvester form. In the SDQS for example: 



C[a] = £ Sy iv[o"] + £ rcst [a] 



where 

^Sylv 

and 

A 00 = 
A n = 



Aa + o-A 1 " 



h °' c 2 

~(-/i/osc cSx 



A)oCoo + o- 00 A 




t 

00 



2x l^j 


•yrnu 




h 






2~ 6 X - 


2h XP 



Ano- u + a u A\ 



ii 



n B 



x 



7mw / 1 

-r[ nB+ 2 lA " 



(3.149) 



(3.150) 



(3.151) 



26 To be precise the preconditioner should be regular only on the image of the Liouvillian. 
27 In this sense the philosophy is: "Invert as much as you can!" 



69 



CHAPTER 3. GENERALIZED MASTER EQUATION 



N = 2 



N = 5 



10 
20 
30 
40 



>i as* 



,20 40 

nz = 246 



N = 10 




100 200 
nz = 3456 



500 
nz = 26006 



1000 



Figure 3.5: Example of matrix representation of the Liouvillean for the TDQS. 
The Hilbert space of the harmonic oscillator has dimension respectively N = 2, 5, 10 
(nz is the number of no-zero elements of the mtrix). The non-zero elements are 
represented by a blue dot. The exponentials of position operator present in the 
tunneling amplitudes make the Liouvillean non-sparse (large squares progressively 
emerging in the figures from left to right). For large N the SVD is not viable and 
we use the Arnoldi iteration scheme to find the null space of the Liouvillean. 



where % is the average occupation number of the energy states of the har- 
monic oscillator in equilibrium with the bath. The rest of the Liouvillean in 
the non-Sylvester form reads: 



r R e*cr u e* - ^(xa 00 p - po~oox) + TI ^ L { c ln B + l)xa m x 
T L e~^a 00 e~^ - ^(xa n p - pa n x) + ^(2n B + l)xa n x 

(3.152) 

We use as the preconditioner for the Arnoldi iteration scheme M. = £gy lv . 
The check of empty intersection between the spectrum of A and —At must 
be done numerically due to the presence of the exponentials of position op- 
erators. 

The stationary solution of the GME represents the starting point for the 
analysis of the properties of shuttle devices that we will present in the fol- 
lowing chapters. The Arnoldi scheme allows the calculation of the stationary 
solution a stat of the GME for a given set of parameters and a fixed dimension 
of the Hilbert space of the oscillator. The number N has been chosen from 
considerations on the phase space distribution for the stationary solution 28 . 



O~00 









28 See next chapter for details on the Wigner distribution and its meaning. 
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Chapter 4 

Wigner function distribution 



The shuttle dynamics has, especially in the single dot device, an appeal- 
ing simple classical interpretation and one can say that the name itself of 
"shuttle" suggests the idea of sequential and periodical loading, mechani- 
cal transport and unloading of electrons between a source and a drain lead. 
Semiclassical models typically treat the mechanical degree of freedom clas- 
sically combined with incoherent sequential tunneling of electrons between 
leads and oscillating dot [1, 2, 3, 4, 5]. Motivated by the small size of the 
oscillations of a nanoscale shuttle, and thus by the possibility to observe 
signatures of quantum dynamics of the mechanical degree of freedom, we de- 
cided, following the suggestion of Armour and MacKinnon [7] , to explore the 
model with a quantized oscillator. Nevertheless we also wanted to keep as 
much as possible the intuitive classical picture and to handle the quantum- 
classical correspondence. The Wigner function distribution seemed to us a 
good answer to all these requirements. It allows a clear visualization of the 
numerical results obtained within the framework of the GME and it shows 
in its equation of motion (the Klein-Kramers equation) an explicit quantum- 
classical correspondence (expansion in powers of H) [34]. 



4.1 A quantum phase-space distribution 

The concept of a probability distribution in phase space is problematic in 
quantum mechanics due to the uncertainty principle. It is not clear for 
example the meaning of the probability that a particle has a well defined 
position and momentum at a certain time and thus it seems difficult to define 
a probability distribution for such "non-observable events". Nevertheless 
one can demonstrate that it is possible to associate with a state defined 
by a density operator p(t) a unique function P(q,p,t) on the phase space 
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(q,p) x [to, oo) that satisfies some desirable properties. We list the required 
properties as stated by Hillary et al. in the review "Distribution functions in 
Physics: fundamentals" [34]: 

1. P(q,p,t) should be real; 

2. P(q,p, t) should reduce to a probability distribution once the conjugate 
variable is integrated: 



dpP(q,p,t) = (q\p(t)\q) 

r (4-1) 

J dqP(q,p,t) = (p\p(t)\p} 
and thus P(q,p) is normalized 

J dp J dqP(q,p,t) = l (4.2) 



3. P(q,p,t) should be Galilei invariant; 

4. P(q,p,t) should be invariant under space and time reflection; 

5. In the force free case P(q,p,t) should obey the classical equation of 
motion: 

§ - -~i <«> 

at m oq 

6. If Pi(q,p,t) and P2(q,p,t) are the distributions corresponding to the 
states Pi(t) and p 2 (t) the following relation should hold: 

dq{q\pm\q) = 2nh J dq J dp P 1 (q,p,t)P 2 (q,p,t) (4.4) 

Condition 6 is a requirement that comes from quantum mechanical consider- 
ations as can be seen from the following consequences. Assume p\ = p 2 = p 
being statistical mixtures of pure states p a with weights w a 

P = ^2 w *Pa (4.5) 

a 

and P(q,p,t) the associated distribution, then from condition 6 we get 
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dq / dp[P(q,p,t)f 



2tt^ 



(4.6) 



This relation is discarding any too peacked distribution. In particular the 
classical distribution P(q,p, t) = S(q — qo(t))5(p — Po(t)) does not fulfill con- 
dition 6. If we assume now that p\ and p2 are pure states associated with 
orthogonal vectors, then: 



dq J dpP 1 (q,p,t)P 2 (q,p,t) = 



(4.7) 



which implies that the two distributions are not positive everywhere. For 
this reason the distribution we are axiomatically constructing will take the 
name of quasi-probability distribution. 

The distribution function that satisfy properties (1-5) or (1-4 and 6) exists 
unique [35], is defined as 



1 f + °° I £ 



p(t) 



and is commonly called Wigner distribution. The parallel with classical 
mechanics can be pushed even further associating to every observable on 
the Hilbert space the corresponding Wigner function: 



A(q,p) 



A 



^+2/ exp It 



Two important relations can be derived from the definitions 



(4.9) 
and (4.9). 



1. The trace of an operator in proportional to the integral on the phase 
space of its corresponding Wigner function: 



1 



2nh 



dq / dpA(q,p) 



(4.10) 



2. Given two operators A and B and the associate Wigner functions 
A(q,p) and B(q,p), the Wigner function associated to the product of 
the two operators F = AB reads [35] 1 : 



F(q,p) = A(q,p) exp (l^J 5 (<? ; 



P) 



(4.11) 



lr rhis relation was first derived by Groenewold in 1946 in the context of an rigorous 
mathematical correspondence between commutator and poisson bracket in the limit h — > 0. 
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where A is the differential operator: 

^ d d d d (4 12) 

dp dq dq dp 

The arrows indicate in which direction the derivative should be taken. A 
better insight in the expression (4.11) is given by the expansion of the expo- 
nential. Up to the first two terms it reads: 



F(q,p) =A(q,p)B(q,p) 

h fdAdB dAdB 
+- ' 



2i \ dp dq dq dp 
+ hV /d^A^B__ 2 ^A_^B_ + d 2 Ad 2 B\ (4.13) 



2i J \ dp 2 dq 2 dpdq dqdp dq 2 dp 2 

h 

+ l 5 

The expansion shows the power of the formalism to treat the quantum- 
classical correspondence: given a length, a mass, and a time scale for the 
system we can rescale the coordinates and the expansion will appear as an 
expansion in h/S sys where S sys is the typical action of the system. Classical 
systems have a large action S sys ^> h and only the first term in the expan- 
sion is relevant. In the opposite limit S sys ~ K the full expansion should be 
considered. We have now the tools to translate into the Wigner formalism 
the expression for the expectation value of an operator: We start from the 
definitions: 

(A(t)) = Ti{p(t)A} = j dqj dpW(q,p,t)ex V (^j A(q,p) (4.14) 

Inserting the expansion in h of the exponential and using the regularity of 
the Wigner distribution at infinity we obtain the "classical" result: 

(A(t)) = J dqj dpW(q,p,t)A(q,p) (4.15) 

A second important consequence of the Groenewold equation (4.11) is 
the quantum-classical correspondence for the dynamics of the distribution 
function W(q,p,t). In the Wigner formalism the equation of Liouville-von 
Neumann reads: 



dW i 



dt h 



(4.16) 
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where H = H(q,p) is the Wigner function associated with the Hamiltonian 
for the system. Due to the fact that the differential operator A is odd with 
respect of the exchange of the two functions on which it is acting only odd 
powers of the expansion in h/S sys survive in (4.16). In the limit h/S sys — > 
we get the classical Liouville equation: 

dW c 

- W = {H,W C } P (4.17) 

where {,}p = — — A are the Poisson brackets and the superscript 

c indicates the classical limit for the Wigner distribution 2 . 



4.2 Klein- Kramers equations for the SDQS 

With the knowledge acquired in the previous section we will now concentrate 
on the translation in the Wigner formalism of the GME for shuttle devices 
that we derived in the previous chapter. We concentrate only on the SDQS 
first because the technique is completely illustrated in this simpler case and 
also because contrary to the single dot device for the TDQS we will use the 
Wigner distribution only as a visualization tool. In the derivation of the 
GME for the single dot device (see §3.4.1) we identified three components of 
the Liouvillean super-operator: 



& -^coh ^driv ^damp (^'1^) 

distinguishing the coherent contribution that describes the evolution of the 
system without the electrical and mechanical baths from a driving contribu- 
tion where the coupling to the leads is inserted and a damping term that 
takes into account the interaction with the dissipative environment of the 
mechanical degree of freedom. We separately treat their translation into the 
Wigner formalism in the following sections. 



4.2.1 Coherent Liouvillean 

We start recalling explicitly the action of the coherent part of the Liouvil- 
lean on the two non-vanishing 3 electrical components of the reduced density 
matrix (<t o and an): 

2 As can be seen from this derivation, in the framework of the Wigner distribution the 
quantum-classical limit attributed to Dirac — ] — ► {,}p for h — > acquires a precise 
mathematical meaning. 

3 For an explanation on the vanishing electrical coherencies ergi and <j\q see §3.3.4. 
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Ceo, ( a °°)=--J m [HoSC, ? 0] i ) (4-19) 

It is easy to verify from the definition (4.9) that for operators which 
are sum of functions of momentum or position operator only, the associated 
Winger functions coincide with the corresponding classical dynamical vari- 
ables. For example the Wigner function of the oscillator Hamiltonian reads: 

H osc (q,p) = ^- + ^mco 2 q 2 (4.20) 

and analogously the electrostatic component becomes eSq. 

As already mentioned the commutator structure of the Liouville-von Neu- 
mann equation implies that only the odd components of the h expansion of 
the Groenewold operator exp(hA/2i) should be retained. Since the third or- 
der in the expansion contains already derivatives in q and p up to the third 
order, and the oscillator Hamiltonian is quadratic in q and p only the first 
order survives and the coherent part of the Liouvilliean acquires the classical 
form 4 : 

^ [Wnj-y m^(q - d)^ - ^ ) ^ 

where the displacement d = is the equilibrium position of the oscillator 
subject to a constant force e£ and we have introduced a Wigner distribution 
for each of the electrical states of the system. In this charge resolved Wigner 
functions one can read the probability density for the quantum dot to have a 
specific position and (average) momentum when in a specific electrical state. 
In this way we can monitor the correlation between charge and position 
(momentum) of the oscillating dot, important to discriminate the shuttling 
regime from other transport regimes. 

4.2.2 Driving Liouvillean 

The driving term of the Liouvillean for the SDQS is a rate equation somehow 
modified to take into account the operator form of the position-dependent 
rates: 

4 This result is more general and resembles the Ehrenfest theorem for Hamiltonian up 
to second order in x and p. In that case one gets classical equations of motion for the 
expectation values (x) and (p) . 
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--driv 



O"00 
0"11 



-^f {e 2x /\ a n } + T L e^ x a 00 e- x / x 



(4.22) 



where {A, B} = + IM is the anticommutator. 

We calculate first the generic form of the symmetric component of the 
driving Liouvillian in the Wigner representation. It is easier in this case to 
start directly from the definition: 



1 



+oo 



-T L , R e^ x I d£(q 

J — oo \ 

--T L , R e^ x W u (q,p,t) 



<Ju{t) 



g+ 2/ 6XP X 



q 



cxp 



h 



(4.23) 



where % — 0, 1 in accordance with the sign of the exponential respectively 
— , + and the subscript L, R in the tunneling rate T. Only the classical 
component contributes to this term since the Wigner representation is just 
the product of the three Wigner functions. 

Quantum correction are present instead in the anticommutator compo- 
nent of the driving Liouvillean. Again we calculate the generic form : 



2nh 



r 



f{e^\o u (t)} 



q + - ) exp ( — 



e T2 "/ A exp 



HA 
~2i 



Wu{q,p,t) + Wu(q,p,t)exp 



h 
HA 



2i 



-^ TVA g(4(f) 2V - teRi) 



2-^ (0^ 



n=0 



{2n)\ \\J dp 



2" pan 



-,2n 



Wa(q,p,t) 



(4.24) 



where because of the symmetry of the anticommutator only even powers 
of the Groenewold operator give non vanishing contributions (line 2 to 3 
in (4.24) ) and they are equal for both terms of the anticommutator. In 
the last line we appreciate the control parameter for the quantum-classical 
correspondence. Especially in the shuttling regime the typical length, mass 
and time scales of the system are the tunneling length A and the harmonic 
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oscillator mass m and inverse frequency Rescaling the phase space 

coordinates to the adimensional X = q/X and P = p/{mujX) we obtain the 
expansion parameter: 



h f Ax 



mcuX 2 \ A 



(4.25) 



that is the ratio between the zero point fluctuation of the oscillator position 
(Ax z ) and the tunneling length (A). The bigger A the more classical is the 
behavior of the harmonic oscillator. 

All together the driving Liouvillean in the Wigner representation takes 
the form: 



(4.26) 



4.2.3 Damping Liouvillean 

The damping Liouvillean is identical for the two electrical components of the 
reduced density matrix and we give a unified derivation of the corresponding 
Wigner formulation. In terms of density matrices it reads: 

£damp(^) = ~[x, {p,a}] - ^(2n B + l)[x, [x,a]] (4.27) 

Once again we apply the Groenewold modified product 5 (4.11) to obtain 
the Wigner formulation. Since the damping Liouvillean is quadratic in the 
operators x and p we expect only the first two terms of the expansion in the 
operator A to contribute. 

C daiap W = ^pW + im huj (n B + (4.28) 

The term of first order in the derivative is called damping term while the 
second diffusion term. One way to understand the meaning of this names is 
to derive this equation for the distribution function starting from the corre- 
sponding Langevin equations describing the brownian motion of a particle in 



5 The product of three operator can be handled easily since one can split the product 
in two parts and the property (A * B) * C — A * (B * C) holds, where we have indicated 
with * the Groenewold modified product. 
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an harmonic potential: 



P 

x = — 
m 



p = —moj 2 x — 7p + £(t) 



(4.29) 



The first term in the RHS of equation (4.28) is related to the velocity depen- 
dent friction (— jp) while the second comes from the stochastic force £(£). 
The analogy with classical brownian motion can be pushed even further con- 
sidering the rescaling: 

P (4.30) 



X 



x 

T 



mcui 



where t = yj -^j{2n B + !)• The damping Liouvillian in terms of the rescaled 
variables reads: 

£d "" W = 7 ^( P + ^) W (431) 

and only the scaling 6 distinguishes the quantum from the classical description 
While in the low temperature limit 2n B <C 1 and the length scale £ tends 
to the oscillator zero point fluctuation Ax z , in the high temperature limit 
h drops out from the equations and t is given by the the thermal length 

\ _ lk B T 

Finally we collect the results of the previous sections and write the equa- 
tion of motion for the SDQS in the Wigner function formulation: 



dW, 



00 



dt 



2 d p d d . ( 1\<9 2 

mu q— - — — + 7— P + imfiuj \ n B + -1 — 



dp mdq dp 



W, 



oo 



l) n (h\ 2n d 2n W ( 



00 



dp 



2n 



d p d 



d 



mu (q - d)- — + 7— p + r/mtuj [ n B + 

dp m dq dp \ 2 / dp A 



\\ d 2 



ii 



^ (2n)\ \\ 

n=0 v ' v 



11 



dp 



2n 



(4.32) 



6 Quantum contribution (e.g regions of negative distribution W(q,p, t) < 0) could be 
contained in the coherence of the initial condition. These will anyway be decohered by 
the damping Liouvillean itself. 
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Figure 4.1: Example of 3-D visualization of the stationary solution of the Klein- 
Kramers equations for the SDQS. The surface represents the quasi-probability dis- 
tribution of the empty dot W^ at (q,p). Coordinate and momentum are expressed 
in adimensional units: X = q/xo and P = pxo/h, where xq = ■sJhjimiS). The 
distribution shows two distinct features. The probability is high near the origin 
of the phase space and represents the oscillator in the ground state. The other 
structure corresponds to the part of the limit cycle drawn by the empty dot in the 
shuttling regime. 
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Investigation tools 



The dynamics of the shuttling devices or at least of the simplified models 
for them that we consider, is all contained in the GMEs for the reduced 
density matrices (3.120) and (3.131) or in the corresponding Klein-Kramers 
equations for the Wigner distributions (e.g. equation (4.32) for the SDQS). 
The coupling to the continuum of states in the electrical and mechanical 
baths introduces irreversibility in the system and its statistical properties 
(represented by the reduced density matrix or Wigner distribution) tend to 
a stationary limit. This does not mean that the system itself is in a static 
condition. In particular an electric current sustained by the bias always flows 
through the device driving the movable dot into a variety of different dynam- 
ical regimes. In order to analyze the properties of these different operating 
conditions of the electronic shuttles we use three investigation tools. We 
study the phase space distribution functions associated with the stationary 
solution of the GME, the corresponding current and the current-noise, taking 
as control parameters the damping rate of the mechanical degree of freedom, 
or the gate voltage of the different dots, or the injection and ejection rates. 
Through the phase-space analysis or the current and noise calculation we gain 
somehow complementary information about the system. These approaches 
are described in detail in the next three sections. 

5.1 Phase space distributions: dots, rings and 
bananas 

The phase space analysis of the shuttle dynamics is based on the Wigner 
representation of the stationary solution of the GME: 




(5.1) 
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where the subscript ij indicates the electrical state that we are considering. 
Since typically we are not interested in coherencies 1 between electronic states 
we set % — j and the resulting Wigner distribution represents the joint proba- 
bility for the system to be in the specific mechanical state (q,p) and electrical 
state i. The total Wigner distribution 

W t T(q,p) = Y,Wt t (q,P) (5.2) 

i 

where the sum is extended to all the electrical states of the system, gives 
information about the mechanical state independently from the electronic 
state of the device. 

5.1.1 Shuttling instability 

For the SDQS we analyze first the behavior of the Wigner distribution as 
a function of the mechanical damping. The mechanical dynamics shows a 
typical feature qualitatively independent from the other parameters of the 
model. At high damping rates the total Wigner distribution is concentrated 
around the origin of the phase space and represents the harmonic oscillator 
in the ground state. While reducing the mechanical damping a ring pro- 
gressively develops and the central "dot" eventually disappears (Fig. 5.1). 
We interpret the ring as a smeared limit cycle trajectory: by reducing the 
damping the equilibrium position becomes unstable and a stable limit cycle 
develops. The threshold for this transition is given by the effective tunneling 
rates of the electrons 2 . 

This interpretation is represented in the following picture: every time an 
electron jumps on the movable grain the latter is subject to the electrostatic 
force e£ that accelerates it towards the right. Energy is pumped into the 
mechanical system and the dot starts to oscillate. If the damping is high 
compared to the tunneling rates the oscillator dissipates this energy into the 
environment before the next tunneling event: on average the dot remains 
in its ground state. On the contrary for very small damping the relaxation 
time of the oscillator is long and multiple "forcing events" happen before the 
relaxation takes place. This continuously drives the oscillator far from equi- 
librium and a stationary state is reached only when the energy pumped per 
cycle into the system is dissipated during the same cycle in the environment. 

It is not difficult to realize that, like for a macroscopic swing, in order 
to sustain the motion one needs coordination between forcing (here related 

1 The probabilistic interpretation of the Wigner distribution would not be acceptable 
for electronic coherencies. 

2 A more precise definition of this electronic time can be found the next chapter. 
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to the electrical dynamics) and oscillations. This coordination is revealed 
by the charge resolved Wigner distributions Woo an d Wu- The ring that 
appears in the total distribution is asymmetrically shared by the empty and 
charged-dot distributions (Fig. 5.1). 



y = 0.25 y = 0.1 y = 0.02 




Figure 5.1: Phase space picture of the tunneling-to- shuttling crossover. The re- 
spective rows show the Wigner distribution functions for the discharged (Wqq), 
charged (W\i) and sum (W to iJ states of the oscillator in the phase space (horizon- 
tal axis-coordinate in units of xq = \Jhj{fujS), vertical axis - momentum in h/xo). 
The values of the parameters are: A = xq, T = 0, d = 0.5xo, P = 0.05u;. The 
values of j are in units ofoo. The Wigner functions are normalized within each 
column. 

The Wigner distributions reveal the charge-position (momentum) corre- 
lation typical of the shuttling regime: for negative displacements and positive 
momentum (i.e. leaving the source lead) the dot is prevalently charged while 
it is empty for positive displacements and negative momentum (coming from 
the drain lead). Self-sustained mechanical oscillations and shuttling electron 
transport are part of the same phenomenon. 
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5.1.2 Classical vs. Quantum Shuttle 

The instability we showed in the previous section was predicted already by 
Gorelik et al. [1] (even if as a function of the external bias) in the form of a 
sharp transition between an equilibrium and a limit cycle dynamics of the os- 
cillator and is called shuttle instability. The quantum mechanical description 
of the harmonic oscillator that we give introduces two effects in the model: 
it smears the transition into a crossover and reveals the possibility to trigger 
the shuttling regime even at zero electric field [8] when classically the elec- 
trical and mechanical dynamics would be decoupled and the device would 
behave as a damped harmonic oscillator that with no feedback influences a 
single electron transistor (Fig. 5.2). 



w 



00 







W 



f i 



IQl 



Figure 5.2: Shuttling at zero electric field. The three Wigner distributions for the 
empty, charged, and sum states are represented in the usual adimensional phase 
space (X,P) (see Fig. 5.1). The parameter values are A = xq, T = 0, d = 0, 
r = 0.05w, 7 = 0. 



The amplitude of the shuttling oscillations of the quantum dot is also 
dependent on the tunneling length and bare injection rate. In particular it 
grows with the tunneling length A and with the inverse of the bare injection 
rate. This behavior can be understood as a self-consistent compensation 




Figure 5.3: Transition from "quantum" to "classical" shuttling. The figure shows 
the total Wigner function for three different sets of parameters. We keep constant 
damping (7 = 0.02a;j temperature (T = 0) and electric field (d = 0.5xoj and 
change the injection rate and tunneling length. From left to right respectively: 
(r = 0.05w, A = xq), (r = 0.05a;, A = 2x ), (T = O.Olw, A = 2x ). 
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to accomplish the non-adiabatic condition that the effective injection rate 
P must be comparable in the shuttling regime with the oscillator frequency 
uo/27T. The larger the amplitude of the oscillation the more classical is the 
behavior of the harmonic oscillator. This concept can be understood also 
from a geometric point of view: the ring gets closer and closer to a circle 
since the width, given by the amount of thermal and shot noise remains 
constant (or decreases as we will see in the coexistence regime semianalytics) 
while the diameter increases. Another natural parameter that reveals this 
quantum to classical transition is the ratio between h and the area enclosed 
by the ring. This tends to zero in the classical limit (Fig. 5.3). 



5.1.3 Different electronic processes in the TDQS 

The Wigner distribution analysis of the TDQS is performed in analogy with 
the SDQS on the states corresponding to the empty and charged central dot. 
Namely we consider: 
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(5.3) 



As suggested by the work of Armour and MacKinnon [7] we consider as 
control parameter the difference between the gate voltages on the left and 
right dot (the so-called device bias AV) . It was found in [7] that the triple dot 
system exhibits different regimes of transport at different device biases. The 
current peaks at AV ~ ntkv were identified as effects of electromechanical 
resonances within the device. Yet, the different peaks may correspond to 
different physical mechanisms - while the peak around AV ~ ftui is mainly 
due to the incoherent oscillator-assisted tunneling the peak at AV « 2f\w 
reveals a clear shuttling component. This finding by Armour and MacKinnon 
based on indirect evidence of parametric dependencies of the current curves 
(e.g. the dependence of the current curve on the tunneling length A) is 
confirmed by the direct inspection of the phase space distribution (Fig. 5.4). 
The half-moon-like shape characteristic for the shuttling transport is only 
present around AV ~ 2u while all other plots show just the fuzzy spot 
indicative of incoherent tunneling. However, our direct criterion for detecting 
the shuttling regime reveals a close similarity between the resonances. For 
increasing injection rate we can see that the shuttling regime gradually sets 
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in also in the vicinity of the first resonance peak. This reveals the cross over 
character of the onset of the shuttling instability found also for the SDQS. 

AV= 0.55 AV= 0.85 AV = 1.3 AV = 1.925 AV = 2.65 
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Figure 5.4: Wigner distributions for the central dot in the charged state. The 
half-moon shuttling signature, always present at the second resonance (AV ~ 2huj) 
progressively emerges also in the first resonance (AV ~ fojj). The other parameters 
are V = 0.757^, A = 5Ax z , x = 7.07lAx z , 7 = 0.0125a;. 



As can be seen in more detail in the work by C. Flindt et al. [33] and in the 
Master thesis by C. Flindt [12] the TDQS exhibits at higher device biases 
another transport mechanism that competes with the shuttling transport. 
For A V ~ 3c<j, 4u inelastic co-tunneling takes place and a clear signature 
of this phenomenon can be seen in the vanishing of the charged central dot 
Wigner function (see Fig. 5 in [33]). 



5.2 Current 

The Wigner distribution functions are a very clean theoretical tool to inves- 
tigate the dynamical properties of shuttle devices, but unfortunately it seems 
very difficult (if not impossible) to access experimentally these functions. For 
this reason, with in mind the distribution function description, we investigate 
the electronic transport properties in the shuttle devices (more realistically 
measurable) starting with the stationary current. 
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5.2.1 Calculation of the stationary current 

We already encountered in section 3.3.5 the definition of right and left lead 
currents: 

I R (t) = N R (t) 

W . W (5.4) 
I L {t) = N L (t) 

where N R (t) (N L (t)) are the total number of electrons (holes) collected in the 
right (left) lead at time t. We have also related these functions to the elements 
of the n(/i)-resolved density matrix and corresponding GME (Eq. (3.54) and 
(3.56)). Since we are interested exclusively in the stationary current 3 we 
analyze now in detail only the electron current in the right lead in presence 
of a shuttle device (SDQS or TDQS). The total number of electrons collected 
in the right lead N R (t) can be written as: 

oo 

N R (t) =Y, nP n{t) (5-5) 

n=0 

where P n (t) is the probability for n electrons to be collected into the right 
lead at time t. We express the probability P n (t) as the trace of the n-resolved 
reduced density matrix over the system degrees of freedom: 

Pn(t) = Tr sys {^(t)} = Tr osc fa*™®} (5.6) 

where % is the electrical state of the device: and in the SDQS % — 0, 1 while 
in the TDQS i = 0, L, C, R. The n-resolved GME (n-GME) comes now into 
play and the current is expressed as a function of (t) : 



oo 

n=0 i 

oo 

= E n E T ^ {A»h[^ (B) (*)]« + AHv W0]i B) + Atamp^W] } 
n=0 i 

(5.7) 

where we have implicitly specified with the position of the square parentheses 
the different components of the density matrix addressed by the different 
Liouvilleans 4 . 

3 The stationary current is equal on both leads due to charge conservation. 
4 Thc damping Liouvillcan acts only in the Hilbert space of the harmonic oscillator and 
does not mix different electronic states or states with different number of particle in the 
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We calculate first the SDQS right (particle) current. The coherent and 
damping part of the Liouvillean have a commutator structure (see Eq. (3.116)) 
and vanish under trace. We are left with the contribution given by the driving 
Liouvillean: 

I R (t) = Tr osc J IV*/* f> [ag-^f) - <7<?(f)] } (5.8) 

I n=0 J 

where we have used the cyclic property of the trace to evaluate the sum 
over the two electrical components. The particular ordering of the position 
dependent injection and ejection operators I^l is irrelevant under trace 
and the sum over the number of electron in the resevoir is formally identical 
to the one obtained for the static single dot (see Eq. (3.54)). We evaluate it 
and obtain: 

I R (t) = Tr osc {T R e 2x / x a n (t)} (5.9) 

With completely analogous arguments one derives the expression for the left 
current: 

I L (t) = Tr osc {V L e- 2x ' x a m {t)} (5.10) 

The stationary current is obtained by substituting in (5.9) or (5.10) the 
stationary solution (j stat of the GME. It is interesting that the stationary 
particle current is written as the average injection (ejection) rate: 

Tr osc {r*e 2 */Vi at } = Tr sys {r i? (x)|l)(l|a stat } = T R 

rstat \< -p 

\ / 

Tr osc {T L e- 2 *l x o^} = Tr sys {r L (x)|0)(0| C T stat } = f L 

(5.11) 

where we have introduced the projectors |1)(1| = qq and |0)(0| = 1 — qq 
on the electronic Hilbert space and we calculate the average injection and 
ejection rates taking into account Coulomb blockade and the single level 
structure. Namely the injection (ejection) rate average is calculated only 
on system states with empty (charged) dot. A special remark concerns cor- 
relations: the shuttling regime is characterized by strong electromechanical 
correlations and they should be recorded in the current since in general: 

J stat = Tr sys {r L (x)|0)(0|<7 stat } + Tr sys {r i (x)a stat }Tr sys {|0)(0| ( T stat } (5.12) 



collector. The coherent Liouvillean is sensitive to the charge state of the system and does 
not care about the particular bath state, while the driving term is mixing electronic states 
and bath states with different number of particles since it describes the electron tunneling 
between system and baths. All of them are local in time (Markov approximation) . 
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5.2.2 Current characteristics of the SDQS 

We study the current as a function of the mechanical damping. The results 
for different values of the bare injection (ejection) rate and of the tunneling 
length are presented in Fig. 5.5. At low damping the current saturates, 
independently from the value of the other parameters, at the "magic" value 
/ pa 0.16a; (numerical approximation of the frequency of the harmonic oscil- 
lator u/2n). Increasing the damping the stationary current drops, more or 
less rapidly, to a plateau, dependent this time on both the bare tunneling 
rate and length. Further increase of the damping does not change the sce- 
nario. The transition between the plateau and the saturation point is faster 
the more "classical" the other parameters are (i.e. small injection rate and 
large tunneling length 5 ). 

Current 

0.2 
0.18 
0.16 
0.14 
0.12 

3 

a> 0.1 

0.08 
0.06 
0.04 
0.02 


0.05 0.1 0.15 0.2 0.25 

y/ co 

Figure 5.5: Particle stationary current of the SDQS plotted as a function of the 
damping rate. The bare electronic rate (T) is given in units of the oscillator fre- 
quency. The tunneling length is given in units of xq = Wh/rnu. The current and 
the damping rate are expressed in adimensional units also by rescaling with the 
oscillator frequency cj. The saturation of the current to the limit of one electron 
per cycle characteristic of the shuttling regime at low damping is visible. The tran- 
sition to the the high damping limit (oscillator in ground state) gets sharper the 
more classical is the dynamics as is indicated also by the parameters (see text for 
explanation). 

If we compare the current results with the Wigner function distribution we 
can recognize a correspondence between the shuttling charge-position (mo- 

5 See previous section for the meaning of the term "classical" associated to a parameter 
configuration. 
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mentum) correlation and the saturation point, as well as a progressive disap- 
pearing of the ring structure in correspondence with the current transition. 
The high damping plateau in the current sets in when the mechanical oscilla- 
tor lays into its ground state and the Wigner distribution function is reduced 
to a fuzzy spot close to the origin of the phase space (Fig. 5.6). 



Current-Wigner 




0.05 M 0.15 0.2 0.25 

| / to 



Figure 5.6: Most "classical" current- damping curve. The Wigner function distri- 
butions in the insets correspond to the three dots indicated on the curve and show 
how the transition looks like in the phase space. 



5.2.3 Electromechanical resonances of the TDQS 

The calculation of the TDQS stationary current follows the same argument 
line we adopted for the SDQS. Also the analytic expression shares essentially 
the same structure, the only difference being the projector that appears in 
the right lead current and the tunneling rate that in this case is not position 
dependent. The electrons tunnel to the resevoir only from the right dot and 
the current reads: 

Jf* = iyTr osc {a**} = r R Tr sys {\R){R\a stat } = T R (5.13) 

Since there are not position operators in the left or right lead current one 
could argue that the stationary current of the TDQS do not record signs 
of charge-position correlation. The solution to this problem comes from the 
structure of the device. Instead of looking at the system-lead tunneling 
current one can concentrate on the tunneling current within the device: in 
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the stationary limit charge conservation ensures they are equal. For example 
the left-center current reads in the stationary limit: 

Itc = Tr osc {it L (x)(a^f - (5.14) 

as can be derived tracing the second equation of (3.131). 

A detailed analysis of the current in the TDQS with the device bias 6 as 
control parameter can be found in the work by Armour and MacKinnon [7] . 
Extension and analysis of the same model in other parameter regimes are 
given in [11] and extensively in [12] and [33]. Qualitatively we can describe 
the current-device bias characteristics as series of peaks (electromechanical 
resonancies) at AV ~ nfru. It is difficult to establish only from the analysis 
of the current the nature of these peaks. Armour and MacKinnon have an 
interpretation of the different mechanisms underlying these resonances based 
on the analysis of the spectrum of the decoupled three dot device + harmonic 
oscillator. This method can visualize very well the basic mechanisms of elas- 
tic or inelastic phonon assisted tunneling. Other information can be gained 
from the behaviour of the current at different resonances as a function of 
the tunneling length or mechanical damping. We acknowledge the validity 
of this analysis but we also think that an independent source of informa- 
tion (such as the phase space or the current-noise) is needed to discriminate 
directly between different operating regimes. In Fig. 5.7 we report a set of 
current curves in which one can recognize the first three electromechanical 
resonances. 

The temperature influences the device characteristics in different ways 
since different competing transport mechanisms (phonon assisted-tunneling, 
co-tunneling, shuttling) depend on temperature in very different ways. A 
generic feature can though be recognized: at finite temperature electrome- 
chanical resonances are found also at negative device biases (see for example 
Fig. 2 in [33]) when electrons can tunnel through the device only by gaining 
energy from the oscillator. This observation suggests the idea of an active 
phonon cooling device 7 useful to prepare for example a macroscopic harmonic 
oscillator close to its quantum ground state more efficiently than by coupling 
it to a "very cold" thermal resevoir. 



6 We recall that the device bias is defined as the difference between the gate voltages of 
the left and right dot and is represented in the model (Eq. (2.9)) as the mismatch between 
the left and right dot energy levels (AV). 

7 This idea was suggested by prof. K. M0lmer during private communications. 
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Figure 5.7: Current as function of device bias in TDQS. Different curves corre- 
spond to different injection rates T. The particle current and the injection rate T 
are scaled with the mechanical frequency uj. The first three electromechanical res- 
onances are visible. The arrows indicate points of minimum or maximum current 
and the corresponding Wigner distributions are reported in Fig. 5.4- The other 
parameters are Vq = 0.757?iu;, A = 5Ax z , xq = 7.071 Ax 2 , 7 = 0. 0125a;. 



5.3 Current noise 



An unequivocal experimental observation of the shuttling transition has not 
yet been achieved. The IV-curve measured in recent experiments on a vi- 
brating C 60 can be interpreted in terms of shuttling [19], but also alternative 
explanations have been promoted [20, 21, 36]. It is therefore natural to look 
for more refined experimental tools than just the average current through 
the device. An obvious candidate is the current noise spectrum [37, 38]. The 
measurement of the noise spectrum or even higher moments (full counting 
statistics) reveals more information about the transport through the device 
than just the mean current. The theoretical studies of the noise have at- 
tracted much attention recently in NEMS in general [39, 40, 41] as well as 
for the shuttling set up [5, 42, 14] in the (semi) classical limit. The quantum 
mechanical theory for shot noise in NEMS that we treat here is also exposed 
in [23] and in a more general form in [33] and is mainly due to T. Novotny 
in collaboration with the author and C. Flindt. 
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5.3.1 The MacDonald formula 

Let us consider the (particle) current flowing from the dot to the right lead of 
the SDQS. This is defined (in the Heisenberg picture) as the time derivative 
of the number operator of the right lead: 



i R (t) = N R {t) = ~[H(t),N R {t)] 

k 

where 

^ = E c U fe , (5.16) 

k 

the symbol • indicates the Heisenberg picture and H is the Hamiltonian 
presented in equation (2.1,2.2). The current-current correlation function of 
right currents reads: 

cwm') = \{{i R ®Mt)}) - (Ut))(i R (t')) (5.17) 

and the brackets ( ) indicate the quantum mechanical average 8 . The noise 
spectrum is defined as the Fourier transform of the correlation function 9 : 

/+oo 
drC lR j R (t + r,t)e itJJT (5.18) 
-oo 

It is possible to demonstrate (see for example [33]) that the zero frequency 
component of the current noise is independent from the particular junction 
considered. We will for this reason speak in general of current noise S(u = 0) 
even if the specific realization of the calculation refers to the right junction 
current. Multitime averages (e.g. (5.17)) can be evaluated in general using 
the Quantum Regression Theorem [27]. Necessary condition is though the 
fact that the operators averaged must belong to the system. This condition 
is obviously violated by the operator /«(£) defined in equation (5.15) that 
contains system and bath operators. An alternative way of calculating the 



8 Since we work in Heisenberg picture the quantum average is given by the expectation 
value on the initial state - i.e. at the time when the Schrodinger and Heisenberg pictures 
coincide. 

9 Some authors define this Fourier with an extra factor of 2 at the numerator. A 
balancing factor of 2 at the denominator is then inserted in the definition of the Fano 
factor. 
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current noise relates it to the probability P n (t) that n electrons have tunnelled 
through the device by the time t. 

We now present this alternative formulation 10 . Let us take the defini- 
tion of noise spectrum (5.18) in the zero frequency limit and consider the 
symmetry relation 



lim C lRtlR {t + r,t)= lim C lR j R {t - r, t) 



t— »oo 



t— >oo 



(5.19) 



and the definition of the current operator (5.15). We can rewrite the noise 
in the form: 



S(0) = lim ({ lim [N R (r) - N R (0)\, ^-[N R (t) - Nr(0)}}) 

t— >oo L r^oo CLt 



d 



(5.20) 



- 2< lim [N r (t) - N R (0)])(-[N R (t) - N R (0)}) 

t^oo CLt 



where we have also used the property: 

lim lim N R (t + r) - N R (t) ~ lim N r (t) - N R (0) 



(5.21) 



where with the symbol ~ we mean that they have the same asymptotic 
behaviour. This property is a direct consequence of a large times stationary 
current. 

Finally, by inserting the definition of transferred charge operator Q R (t) = 
N R (t) - N R (0) into (5.20), we obtain: 



5(0) 



lim — 

t— »oo at 



(Ql(t)) - (QrW 



(5.22) 



In the basis in which the transferred charge operator is diagonal we can ex- 
press the quantum expectation values 11 of Q and Q 2 using the corresponding 
diagonal density matrix 12 P n (t) and we obtain: 



S(0) = hrn^ j f [ J2 n 2 P n (t) - ( nPnit)) 2 } (5-23) 



n=0 



ra=0 



known as the MacDonald formula [43] . This formula represents the starting 
point for the calculation of the noise in the shuttle devices. We analyze in 



10 One of the main reasons for the extensive derivation of the n-GME that we gave in 
chapter 3 is this formulation of the SDQS current noise. 

11 In what follows we assume the equivalence between the Heisenberg and Schrodingcr 
picture and evaluate the averages in the latter representation. 

12 Note that in this basis the density matrix is diagonal since the state of the right bath 
is a statistical mixture of states with different particle numbers. 
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detail the SDQS case. The TDQS case is the main subject of a detailed article 
[33] and of the Master thesis of C. Flindt [12] and we direct the interested 
reader to that literature. 



5.3.2 Current noise in the SDQS 

The shot noise is completely determined by the evaluation of large times 
asymptotics of the functions: 



A(t) = Y,nP n (t), B(t) = J2nPn(t), C(t) = J> 2 P n (t) (5.24) 

n=0 n=0 n=0 

that appear in the MacDonald formula if one explicitly calculates the time 
derivative. We recognize in A(t) the average number of electrons N R (t) 
tunnelled to the right lead by the time t and in B(t) the right lead current 
Iii(t). From the previous section we know the asymptotic behaviour of the 
the current at large times: 

Km Jfl(t) = Tr osc {T R e 2x ' x aff} (5.25) 

t^oo 

Using the n-GME we can express also the C function in terms of the reduced 
density matrices 



c(t) = x> 2 p n (t) = xy^Tr^/vir^ - 

(5.26) 

2^Tr osc [r^Mf] + r*Tr osc [e b V] 



n=0 n=0 

oo 



n=0 



where we have used the boundary condition crj^ = due to the high bias 
limit that does not allow electrons to enter from the right lead. It is con- 
venient for the calculation of the t — > oo limit in (5.23) to introduce the 
generating functions 13 : 

oo 

F«(f;*0 = £<72°(*)* B (5-27) 

n=0 



13 These are in fact an operator- valued generalization of the generating function concept. 
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with the properties: 

Fu(t; 1) = <Ju(t) 



9 P/ ,, _^__(»),^ (5-28) 



n=0 



The equations of motion for F u (t;z) are easily derived from the GME for 
the n resolved density matrix (3.98) by multiplying both sides with z n and 
summing over n. The resulting Liouvillean for the F's is very similar to the 
one we calculated for the density matrix 14 : 

d i 

O^F 00 (t; z) = --[H osc , F 00 (t; z)\ + £dam P -^oo(^ z ) 

- -j-{e~^,F m {t- z)} + zT R e^F u {t; z)e* 



(5.29) 



= £ 00 F o(t; z) + zC 01 F n (t; z) , 

i 
— F n (t] z) = --[Hose - eSx, F n (t; z)] + C damp F u (t; z) 

- y{ex F u (f; z)} + T L e~* F 00 (t; z)e~x 
= C w F 00 (t; z) + £uF n (t] z) , 

where we have introduced the block structure of the Liouvillean (super)operator 
£ = ). Using the F's and (5.26) the shot noise formula (5.23) can be 

rewritten as 



5(0) 



(tt osc [r R eX (^F u (t; z)\^ + F n (t; 1) 



-2Tr n 



2x r d x — > 

T R e-F u (t; 1)J x Tr osc [— z) 



The Laplace transform of (5.29) yields 



(Foo 


>, 


z^ 




fs - Coo ~ 




>, 




1- 


I -^io s 



i=0 



-01 

Cn 



2=1J 



(5.30) 



t— »oo 



(5.31) 



where flf lt (z) = J2 n a^\o)z n depends on the initial conditions. Defining the 
resolvent Q(s) — (s — £) _1 of the full Liouvillean we arrive at the following 

14 The only small (but important) modification is the z in the driving term of the first 
equation. It corresponds to the change of the electron counting exponent to n — 1 in the 
original n-GME (Eq. (3.98)). 
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expressions for the operator-valued generating functions F^s 

& 11) = sis) (S (5 - 32) 



and their derivatives 15 d z Fu 



evil ^^)f:SV^)^ lt(i; 



d_(F 00 (s;z) 

dz \F u (s;z)J z _ x y ™ {0 ; V<V ' * V "' V/n nit (l), 

(5.33) 

In order to extract the large-t behavior we study the asymptotics of the above 
expressions as s — > 0+ 16 . This is entirely determined by the resolvent C/(s) 
in the small- s limit. Since £ is singular (recall Ca stat = 0) the resolvent is 
singular at s = 0. To extract the singular behavior we introduce the projector 

/ _stat \ 

V on the null space of the Liouvillean: V* = ( ^at ) Tr sys («). We also need 

the complement Q = 1 — V . The projectors V, Q and the Liouvillean C 
fulfill the relations: 

<PjC = J C<p = 

5.34 

C = QCQ K } 

These relations have a simple geometric interpretation. The Liouvillean (su- 
per)operator annihilates the null vector component of the operator on which 
it is applied. This component is instead the only one preserved by V . Every 
vector is thus sent to zero by the combined action of C and V in whatever 
order. The second relation in (5.34) reflects the complementary behaviour of 
the projector Q that extracts the operator components which are also pre- 
served by the Liouvillean and does not affects the action of the latter. The 
resolvent can be expressed as: 

Q(s) = (sV + sQ-QCQ)- 1 

111, (5.35) 
= -V + Q -Q^-V-QC-'Q 

s s — L s 

in leading order for small s. The object QC~ l Q (the pseudoinverse of C) is 
regular as the "inverse" is performed on the Liouville subspace spanned by 
Q where £ is regular (no null vectors). 



15 It can be useful in the calculation to remember the matrix identity : J| [A 1 (z)B(z)] = 
A-\z)f z B(z)-A-^l z A(z)[A-'B(z)}. 

16 Given a function f(t) and its Laplace transform f(s) the following identity holds 



lim t ^ +00 f(t) = lim s ^ + sf(s). 
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We can now calculate the s — > 0+ limit of (5.32) and (5.33) and via 
inverse Laplace transform find their large t asymptotics: 



Foo(t;l) 



t—>co 



a 00 

init 
°11 



a 00 

stat 
°11 



(5.36) 



and 

d_ (F 00 (t;z) 
dz \F n {t;z) 



2=l,t— ►«) 



°oo 

^.stat 
°11 



(t/ stat + C init ) - QC^Q 







(5.37) 

where C imt depends on initial conditions. We insert these in (5.30) and 
obtain: 



S(0) = 2Tr osc 
- 2Tr OJ 
"I - Tr osc 
-2Tr os 

+ 2Tr os 



T R e 2x/x alf\tI stat + C 



ini^ 



t— >oo 



T R e 2x / x 
T R e 2x/x aff 



Y R e 2x/x af^ 



QC~ l Q 



fV R e x l x a 



stat 
11 



11 



1 

x Tr osc [^* at (^ stat + C init ) 



i=0 

i 



T R e 2x / x aff xTr osc ^ 



i=0 



t^oo 

T R e x / x al\ at \ 

1 



(5.38) 



The first term cancels identically the fourth and the noise (as it should) is 
independent from the initial conditions. Also the last term vanishes since 
the linear form Tr osc ^ = Tr sys tries in vain to extract the component in the 
direction of null vector for the Liouvillian which is absent after the projector 
Q has been applied. 



5.3.3 The Fano factor 

We choose to represent the noise of the SDQS by the adimensional quantity 
F = S(0)/I st£lt called Fano factor. In terms of the stationary solution of the 
GME the Fano factor for the SDQS reads 



F = 1 



Jstat llosc 



2x 

e * 



ii 



(5.39) 
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as can be proved dividing (5.38) by the stationary current J stat . For the 
numerical evaluation of the Fano factor we introduce the auxiliary quantity 

S = Q£- 1 Q^ ef < tef ) (5.40) 
which satisfies the equation. 

/ V ^x/X^stai, „x/A^_stat rstat ^.stat \ 

CX=F" " u t^T 1 °" (5.41) 



Once the stationary solution of the GME <r stat is known, equation (5.41) is 
just an inhomogeneous (linear) equation for E. Once again the dimension of 
the matrix representation of the Liouvillean (super)operator C makes the nu- 
merical problem non-trivial. The solution comes from a reformulation of the 
ideas and concepts of the Arnoldi iteration scheme applied to inhomogeneous 
equations called Generalized Residual Method (GMRes) [31]. The problem is 
to find, for a given vector b in the Liouville space, the solution x of the equa- 
tion Lx = b. We start from a Krylov space JCj(L, r ), where r = Lx — b 
and x represents the initial guess for the solution. GMRes finds the vector 
Xj G /C,-(L,ro) that minimizes the norm of the residual r, = b — Lx.,. We 
represent the vector x^ in terms of j coordinates of the Krylov space and 
the matrix of the basis vectors: = x + QjVj. We shift the minimum 
problem from the Krylov space /C- to 0: 



min ||b — Lx||2 = ||b — Lxj||2 



= ||b-L(x + Q j v j )|| 2 
= || r - LQ^Vjl^ 

= || r - Q, + iH,v,|| 2 = min ||r - Q, + iH,-v| 



(5.42) 



where we have used the fundamental relation (3.142) to get rid of the Li- 
ouvillian operator. Since the Krylov space is constructed from r the first 
column of Qj+i is jr^V and we thus have: 

ll r o - Qj+iHjVj|| 2 = ||Qj+i(eir - Hjv^y = ||eir - Hj-VjI^ (5.43) 

with r = ||r || 2 and d = [1,0,..., 0] T G C (j+1)xl . The problem of mini- 
mizing ||eir — HjVj|| 2 can be solved using Q-R-decomposition of Hj. This 
method is fast since the dimension of Hj is of the order of the Krylov space di- 
mension. The method is, like Arnoldi scheme iterative and the convergence is 
strongly dependent from the problem at hand. The preconditioning is again 
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a good method to cure non-convergent problems. We used in this case the 
same Sylvester-like preconditioner that we introduced for the calculation of 
the stationary solution of the GME in Sec. 3.5. 

Before presenting the results of the numerical evaluation of the Fano 
factor in the SDQS we want to discuss its typical behaviour for simpler 
systems. This will provide us with some expectations and ideas to interpret 
the numerical results. 

For a system that can be treated in the Landauer-Biittiker formalism 
(see for example [44]) the scattering states and the transmission probabilities 
between in and out states define the transport properties. In particular the 
average (particle) current through the device is given as: 

, 2e, N 



h 



-Vj^Tn (5-44) 
=i 

where T n is the transmission probability for the nth conducting channel and V 
is the bias across the system. The zero frequency current-current correlation 
function reads: 

2 N 

S(u = 0) = -£v ^T„(l - T n ). (5.45) 

71=1 

For these systems the Fano factor takes the form: 

f= m = tz.m-t.) ( , 46) 

which is a number between and 1. For small transition probabilities (T n <C 
1) the Fano factor tends to 1. The tunneling event is so rare that is reasonable 
to think that there is no correlation between two subsequent tunneling events. 
The Fano factor 1 is also called Poissonian since it is possible to demonstrate 
that in that regime the number of electrons N(At) transmitted in a time 
interval At is distributed according to a Poisson distribution: 

P{N(At) = k} = ^IV AA * (5.47) 

where A is the intensity of the process. Davies et al. [45] described the trans- 
port process for a single channel as a classical stochastic process of particles 
with a definite probability T of going through a barrier. It is a classical point 
of view that condenses the dynamics into a probability distribution of success 
for a tunneling event. In their analysis the Fano factor naturally appears as 
a measure of the randomness of the process. Totally uncorrelated tunneling 
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events give a Fano factor 1. Resonant tunneling through a symmetric double 
barrier give Fano factor F — 1/2 whereas a Fano factor F = corresponds 
to transmission with no randomness. 

The strong position-charge correlation of the shuttling regime suggests 
that the tunneling events are almost deterministically determined by the 
mechanical dynamics. Consequently we expect that the low degree of ran- 
domness in the electron transport is reflected in low Fano factors. This con- 
jecture was recently confirmed by Pistolesi [5] that has predicted a vanishing 
Fano factor for a driven classical charge shuttle at large amplitudes. 

Landauer-Biittiker formalism can treat only sub-poissonian noise since 
F < 1 in (5.46). However, in interacting systems, general theoretical pre- 
diction and numerical calculation have demonstrated the existence of super- 
poissonian noise. The Fano factor may even diverge as discussed by Blanter 
and Biittiker [37]. 

In figure 5.8 we present the Fano factor as a function of the mechanical 
damping 7 for different values of the bare injection rate Y and tunneling 
length A. We recognize common features in the three curves. At high damp- 
ing the Fano factor is of order 1 and (at least for the "most classical" set 
of parameters A = 2xq and V = 0.05cj, 0.01a;) close to what we expect from 
a resonant tunneling system. We will see that the discrepancies from the 
symmetric double barrier are due to the quantum fuzziness in the position 
of the dot that influences the injection and ejection rates and to the charge 
dependent equilibrium position of the dot. Diminishing the damping the 
Fano factors encounter a more or less pronounced maximum and drop fi- 
nally to very low values for small damping. The maximum at intermediate 
damping rates is more pronounced and sharper the more classical are the 
parameters and can reach values of F xs 600 for the most classical case. A 
comparison with the current curves (fig. 5.5) shows that the peak in the 
Fano factor corresponds to the transition region from the tunneling to the 
shuttling current. Similarly the correspondence can be established also with 
the Wigner function distribution: the region of damping in which tunneling 
(dot in W tot ) and shuttling (ring in W tot ) features coexist is associated with 
a super-poissonian Fano factor. The very low (F ps 0.01) Fano factors for 
low damping are a signature of the deterministic transport that takes place 
in the shuttling regime. It is interesting to note that this regularity persists 
also deep in the quantum regime as can be seen for T = 0.05u; and A = xq. 
The relative uncertainty in the amplitude of the oscillation (see Fig. 5.1) 
does not seem to influence the current noise 17 . 

17 It must be noticed that this behaviour is fragile and strongly dependent on the degree 
of unharmonicity of the oscillator. The current is given by the mechanical oscillation 
period that, for a harmonic oscillation is independent of the amplitude. Thus amplitude 



101 



CHAPTER 5. INVESTIGATION TOOLS 




Figure 5.8: Fano factor for the SDQS vs. damping 7. The mechanical dissipation 
rate 7 and the electrical rate T = Tl = are given in units of the mechanical 
frequency uj. The tunneling length A in terms of xq = yjh/ (raw), The other 
parameters are d = 0.5xo and T = 0. The very low noise in the shuttling (low 
damping) regime is a sign of ordered transport. The huge super-poissonian Fano 
factors corresponds to the onset of the coexistence regime. 



fluctuations do not influence the electrical dynamics. Increase of the Fano factor in the 
shuttling regime has been reported for unharmonic potentials [42] . 
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The three regimes 



The shuttle devices are characterized by a strong interplay between electrical 
and mechanical degrees of freedom, the electrostatic force modifying the me- 
chanical dynamics of the oscillating dot and the position dependent tunneling 
amplitudes representing the back-action of the mechanical degree of freedom 
on the electrical dynamics. Despite the complexity expected because of the 
non-linear couplings, the dynamics of shuttle devices (at least in the SDQS 
realization) can be classified into three operating regimes defined by different 
relations between the typical time scales in the device. We make in the next 
sections a qualitative analysis of these regimes leaving to the next chapter the 
definition of the simplified models and the quantitative comparison between 
full and approximate description in terms of the investigation tools that we 
introduced in the previous chapter. 



6.1 Tunneling 

The tunneling regime is the high damping regime in which the relaxation 
time of the mechanical oscillation (inverse of the damping rate 7) is much 
shorter than the typical electronic time (inverse of the average tunneling rate 
f). We also assume that the oscillator dynamics is under-damped: 

u » 7 » f = Tr sys {T R (x)\l)(l\a stat } (6.1) 

The dynamics of the device can be described by a four steps cycle (see Fig. 
6.1): ' 

1. The dot is empty and in its mechanical ground state. We represent it 
in the phase space at rest in the origin and neglect in this simplified 
description the quantum indetermination and/or thermal noise that 
would blur the mechanical phase space distribution. 
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Figure 6.1: Schematic picture of the tunneling regime. In the scheme we represent 
the motion of quantum dot in the phase space of the oscillator (position X and 
momentum P) and at the same time we keep track of the charge degree of freedom. 
The four steps on the cycle are described in detail in the text. 



2. A tunneling event (from the source lead) and consequent electrostatic 
force e£ perturbs the mechanical state of the dot which gets a positive 
position and momentum and starts to oscillate due to the harmonic 
restoring forces. The amplitude of the dot oscillations diminishes due 
to the dissipative environment. The 1) — > 2) transition takes place 
randomly but with a defined rate T L : the average injection rate 1 . 

3. The dot is at rest in the charged state equilibrium position (d = - 2 t-) 
given by the combination of the harmonic potential and the electro- 
static force. The transition 2) — > 3) takes place in the short relaxation 
time I/7. 

4. A second tunnelling event (now towards the drain lead) brings back the 

1 In the choice of as injection rate implies that we are taking into account the average 
occupation of the dot. In other words we are treating the stationary state. See the next 
chapter for more details. 
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dot into the original harmonic potential. The empty dot (distribution 
function) spirals towards the origin of the phase space and the cycle is 
closed. The time scale of the 3) — > 4) transition is given by the inverse 
ejection rate 1/Tr while the mechanical relaxation 4) — > 1) takes place 
within a time 1/7. 

Since the mechanical damping rate is much larger than the average injection 
(ejection) rate the system will stay most of the time in the configuration 1 
and 3. A description with a coarse grained time do not see the mechanical 
transients 2 and 3 and is given in terms of a two state model: empty dot in 
the origin and charged dot in the shifted equilibrium position. 

6.2 Shuttling 

The shuttling regime is characterized by a non-adiabatic interplay between 
electrical and mechanical degrees of freedom of the device. In particular the 
average injection and ejection rate exactly matches the mechanical frequency 
of the oscillator. The mechanical relaxation rate is much smaller than the 
mechanical frequency: 

7 « £ = r (6.2) 

In the shuttling regime the quantum dot oscillates in a self-sustained stable 
limit cycle: all the energy pumped per cycle into the system by the electro- 
static force is dissipated in the environment during the oscillation itself. Also 
the shuttling regime can be described by a four step cycle (see Fig. 6.2): 

1. The quantum dot, close to the source lead, is charged with an extra 
electron. The Coulomb blockade prevents any other charging event. 
The high chemical potential on the left lead makes the charging inco- 
herent and irreversible. 

2. Under the combined effect of the mechanical restoring force and the 
electrostatic force the dot is accelerated towards the right lead. The 
tunneling parameters have been chosen such that the dot is effectively 
decoupled from the left and right lead while close to the center of the 
oscillation. Thus the electrical state of the device remains unchanged 
along most of the oscillation trajectory. 

3. The dot is close to the right lead and the very high ejection rate makes 
the electron unloading practically deterministic. 
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4. The empty dot returns towards the left lead under the action of the 
restoring force of the harmonic oscillator and completes the cycle. 

This sequential "classical" description of the shuttling regime is motivated by 
the phase space analysis of the previous chapter. The current saturation of 
one electron per cycle and the low noise seem to confirm this interpretation. 
At least for the most classical parameters range we will see that the dynamics 
is well described by a "trajectory" in the phase space of the system where a 
mean field charge coordinate has been added to the mechanical position and 
momentum 2 . 



1) 




2) 



r= 



CD 



271 



V 



X 



X 





P 



4) 



3) 



Figure 6.2: Schematic picture of the shuttling regime. 



2 See section 7.2 
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6.3 Coexistence 

At intermediate damping rates the shuttle device shows bistability. Under 
these conditions neither the tunneling or the shuttling are globally stable 
regimes and the noise present in the system (shot noise of the charge trans- 
fer, thermal noise, quantum mechanical indetermination) causes the random 
switching between the two processes. The bistability is represented in the 
stationary distribution function by the coexisting tunneling and shuttling fea- 
tures. We showed that tunneling and shuttling regimes are associated with 
different average currents. It is not difficult to imagine that the current of 
the coexistence regime is an average of the shuttling and tunneling currents 
weighted on the probability for the device to be in one or in the other regime. 
The slow switching between these two currents generates the peculiar huge 
super-poissonian Fano factor that also characterized the coexistence regime. 
It is possible to demonstrate that the zero frequency current noise associated 
to a slow dichotomous switching between two states with average currents 
I + and J_ reads: 



where a (/?) is the switching rate from state 1 to 2 (2 to 1). It is impor- 
tant that the average transition time between shuttling and tunneling is the 
longest time scale in the system dynamics. This fact ensures: 

1. Many tunneling or shuttling events before the regime transition take 
place and thus there exists a well defined average current for the two 

"states"; 

2. A huge Fano factor F pa I sh /a. 

The time resolved dynamics for the coexistence regime is depicted in 
figure 6.3 where the transition rates between the two regime are represented 
respectively by Y out (tunneling — > shuttling) and Y in (shuttling — > tunneling). 
The names {in and out) associated to the transition rates are motivated by 
the dynamics of the mechanical amplitude during the transition. The average 
oscillation amplitude are very different in the shuttling and tunneling regime 
and can be used to identify the two states. In fact we will describe the 
switching dynamics in terms of a bistable potential with the amplitude as 
the effective coordinate 3 . 



3 See section 7.3 for details. 



5(0) 



a(5 



(6.3) 
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Figure 6.3: Schematic picture of the coexistence regime. 
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Simplified models 



We qualitatively described in the previous chapter three possible operating 
regimes for shuttle-devices. They represent for the SDQS the whole scenario 
of possible dynamics. Part of the complexity of the TDQS is also captured in 
this scheme 1 . The specific separation of time scales allows us to identify the 
relevant variables and describe each regime by a specific simplified model. 
Models for the tunneling, shuttling and coexistence regime are analyzed sep- 
arately in the three following sections. We also give a comparison with the 
full description in terms of Wigner distributions, current and current-noise 
to prove that the models, at least in the limits set by the chosen investigation 
tools, capture the relevant dynamics. 

7.1 Renormalized resonant tunneling 

The electrical dynamics is definitely the slowest in the tunneling regime since 
the mechanical relaxation time (much longer in itself than the oscillation pe- 
riod) is much shorter than the average injection or ejection time. We already 
noticed (sec. 6.1) that, because of this time-scale separation, the observation 
of the device dynamics would most of the time show two mechanically frozen 
states: 

0. Empty dot at rest in the oscillator equilibrium position 

1. Charged dot in the shifted equilibrium position of the harmonic oscillator 

perturbed by the constant electrostatic force e£. 

We combine this observation with a quantum description of the mechan- 
ical oscillator and possible thermal noise under the assumption that the n- 
resolved reduced density matrix of the device can be written in the form: 

1 For a further insight into this particular model see [7, 12, 33]. 
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a ( n\t)=pi\t)a th (e£) 

where 

e -f)(H osc -Fx) 

a ^ = Tr osc \pm^=m] (7 - 2) 

is the thermal density matrix of a harmonic oscillator subject to an external 
force T . The functions p^(t) and Pu\t) represent the probability to find 
the system respectively in the state or 1 that we described above when 
n electrons have tunnelled through the dot. The sum over all possible de- 
vice and bath configurations is Yin Poo W + Pii (0 = 1 & t all times. The 
equations of motion for the probabilities p^ (t) can be derived by inserting 
the assumption (7.1) in the n-GME (3.98) and taking the trace over the 
mechanical degrees of freedom 2 : 

. (n) fi (n— 1) fi (n) 

Poo = r RPn -^lPoo 

.(n) f, (n) f, (n) v'^' 

where 

f L = Tr osc [r L e-x <Tth (0) 

~ r 2a; 

T/j = Tr osc T R e~a th (eS) 



(7.4) 



are the injection and ejection rates renormalized to take into account the 
quantum/thermal distribution of the harmonic oscillator position. The mas- 
ter equation (7.3) summed over the number n of electrons collected in the 
right lead reads: 



Poo = r R p u - T L p 00 

fin = TlPoo - fflPn 

where we have introduced the occupation probabilities of the state % pa = 
Yn°=oPif' Equations (7.5) describe the dynamics of a resonant tunneling 
device. The effects of the movable grain are contained in the effective rates 



2 The damping Liouvillean £<jamp is missing in 3.98. Anyway, together with the coherent 

ltribr 

under the trace Tr osc because of their commutator structure. 



(n) 

Liouvillean £ co h, it does not contribute to the equation of motion for the p u . They vanish 
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7.1.1 Electrical rates 

The renormalized electrical rates can be calculated analytically as functions 
of the model parameters. The calculation gives the exact contribution of the 
quantum mechanical and thermal noise to the electrical dynamics. We first 
consider the limit of very low temperatures k R T <C frw. Only the ground 
state of the harmonic oscillator is occupied and the renormalized injection 
rate reads: 



f L = Tr osc [|0)(0|r L e- 2 ^ A ] = r L (0|e^|0> 



d£ exp 



m 2 



(7.6) 



where |0) is the vector representing the ground state of the harmonic oscil- 
lator and x = yjh/{rnw). For the calculation of the renormalized ejection 
rate it is useful to introduce the displacement unitary operator: 



with the properties : 



D(l) = (7.7) 
D(l)dD\l) = d-l^£ 

D(l)$D\l) = d) - lyj^ ( 7 - 8 ) 

D(l)xD\l) = x — l 

where x,p are the position and momentum operators and d',d are the cre- 
ation and annihilation operators for the harmonic oscillator. In terms of the 
displacement operator D the density matrix Oth^) can be written: 

M^) = D (£,) a th (0)^t (_£,) (7.9) 
Only the displaced ground state contributes to the renormalized ejection rate: 



f R = Tr osc [D(d)\0}(0\D\d) T R e- 2 ^ x ] 

= r fl (0|D t (d)ex£)(d)|0) (7.10) 

„ 2d , , 2x , . „ 2d ,( x 0\ 2 

= r R e-(0|e-|0) = T R e- + y-) 

3 These properties can be easily derived from the definition of the operator D(l) and 
the Baker-Hausdorff-Campbell formula [25]. 
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where d = 

As expected the ejection rate is modified by the "classical" shift d of the 
equilibrium position due to the electrostatic force on the charged dot. Some- 
how unexpectedly both rates are also enhanced by the quantum uncertainty 
in the position present in the oscillator ground state. The relevance of the 
quantum correction is given by the ratio between the zero point position dis- 
persion xo and the tunneling length A: the smaller the ratio the smaller the 
correction. 

The calculation of the finite temperature rates is problematic due to the 
presence in the trace of the product of the exponentials of two non commu- 
tative operators (H osc ,x). We solve the problem using the low temperature 
calculation just performed and some symmetry arguments. First we notice 
that the density matrix a th is the stationary solution of the GME 4 : 



* = --[H osc ,a] - -lx,{p,*}] - — 



/ s 1 



[x,[x,a\\ (7.11) 



We express the trace in the position (generalized) basis \q) and the renor- 
malized rates take the form: 

_2q 

dqP t h{q)e * 

(7.12) 

2d I , . 2q 

Y R = T R e~ I da P.u(n)e~ 



J dqP th (q)e £ 



where we have introduced the coordinate probability distribution Pth(q) = 
(q\crth\q)- We encountered this distribution already in section (4.1) in the 
axiomatical definition of the Wigner function. We required specifically: 



/ 



d P W(q,p,t) = (q\a\q) (7.13) 



The problem is thus shifted to the calculation of the Wigner function asso- 
ciated with the thermal density matrix a t h- The dynamics of the Wigner 
distribution is given by the Fokker-Planck equation 5 related to the GME 
(7.11): 



dW 



x— - p— + 1— (p + — 

dP dX uodP\ dP 



W (7.14) 



4 We have assumed this property from the very beginning. For a detailed proof of this 
statement sec for example 11 The Theory of Open Quantum Systems" by H.-P. Breuer and 
F. Petruccione [46]. 

5 The Klein-Kramers equations of second order are also called Fokker-Plank equations. 
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where we have defined the adimensional variables: 



X = T p = £i <" 5 > 



and £ = y ^b^-B + -*-)■ ^he Fokker-Planck equation (7.14) does not de- 
pend on the temperature. The stationary solution is the Gaussian 6 in X and 
P 

1 / X 2 4- P 2 \ 

W S ^(X, P) = — exp ^ ±— J (7.16) 

with no parameters. The position distribution function P t h(<z) thus depends 
on the temperature only in the scaling factor and has the same functional 
form (Gaussian) of the zero temperature case: 

Finally we insert (7.17) into (7.12) and perform the gaussian integrals to 
obtain the finite temperature rates: 



f L = r L e 2 (^) 2 

(7-18) 

In the low temperature limit equations (7.18) obviously reduce to (7.6) 
and (7.10). The high temperature limit k R T 3> ho is also interesting: 



f L = r L e 2 ( A *) 

(7.19) 

Y R = Y R e" +2 y * ) 

In this case the renormalization is the more effective the larger is the thermal 
length A t h = \J k B T / {muj 2 ) compared to the tunneling length A. 



6 This statement can be easily verified by substitution. In order to derive it one con- 
siders that the Gaussian in P is the solution for the dissipative part of the Liouvillean 
(proportional to 7) while the coherent part is solved by spherically symmetric distributions 
since in polar coordinates (A, <j>) we have Xdp — Pdx — > — 9$. 
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7.1.2 Phase-space distribution 

The phase space distribution for the stationary state of the simplified model 
for the tunneling regime is built on the Wigner representation of the thermal 
density matrix <x t h and the stationary solution of the system (7.5) for the 
occupation p u of the electromechanical states i. The latter is easily expressed 
in terms of the renormalized rates T L R : 



stat _ 

Poo — ^ 



R 



stat _ 
Pll — ^ 



T L + T R 



(7.20) 



R 



We have already met the Wigner distribution corresponding to the thermal 
density matrix <7 t h in Eq. (7.16). We now rewrite it in terms of the canonical 
variables (q,p): 



<J th (0)-^W th (q,p) = 



2-KmojP 



cxp 



-) 2 + 

U 



v y 



(7.21) 



The calculation of the Wigner function for the displaced thermal distribution 
Cth(e£) can be performed using the properties of the displacement operators 
D: 



o~th(e£) 



2irh 
1 



D(d)a th (0)D\d) 



q + ^ ) exp I -p£ 



fi 



= W th (q-d,p) 



d ~2 



^th(O) 



9 ~ d + o ) ex P I jP^ 



(7.22) 



where we have used the property D{l)\a) = \a — I) that we can be proven as 
follows: given \a) eigenvector for the position operator x with eigenvalue a 



xD(l)\a) = D(l) D\l)xD(l)\a) 

= D(l)(x-l)\a) = (a - l)D(l)\a) 



(7.23) 



The stationary solution of the Klein-Kramers equation (4.32) in the tun- 
neling regime limit thus reads: 



TT^stat 
vv 00 



(q,p) 



WH at (q,p) 



R 



=rW th (q :P ) 



(7.24) 
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Figure 7.1: Wigner distribution functions for the SDQS in the tunneling regime, 
calculated numerically using the full description. Coordinates (horizontal axis) 
are measured in terms of xq, momenta (vertical axis) in terms of mujXQ. The 
crosses indicates the cuts in the phase space along which we perform the comparison 
between numerical and analytical solution shown in figure 7.2. The parameters are: 
T L = T R = 0.01a;, 7 = 0.25a;, d = 0.5x , A = 2x , T = 0. 




Figure 7.2: Comparison between the numerical and the analytical results for the 
Wigner distribution functions. The squares (circles) are numerical results for the 
parameters mentioned in figure 7. 1 in the charged ( empty ) dot configuration, and 
the full lines represent the analytical calculations. We also plotted with dots the 
numerical results for T = 0.001a; . 

The stationary distribution of the tunneling model is completely deter- 
mined by the length I and associated momentum mui, the equilibrium po- 
sition shift d and the tunneling length A. The electronic and mechanical 
relaxation rates 7 and T^r drop out from the solution and only set the 
range of applicability of the simplified model. 

In figures 7.1, 7.2 and 7.3 we compare the Wigner functions calculated 
both analytically and numerically in the tunneling regime. They show in 
general a very good agreement. In particular we notice in Fig. 7.2 that the 
convergence to the tunneling regime is not yet complete for 7 = 0.25a; and 
^l,r — 0.01a;. To check the theory we diminished the bare electrical rate fur- 
ther (the quantum corrections are negligible since A = 2x , the classical shift 
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can be taken into account and gives a factor e ~ 2.7 which does not change 
the analysis) and kept the mechanical damping fixed to "maintain" the con- 
dition 7 C w. For T = 0.001c<j we got perfect agreement. We also analyze 
the temperature dependence of the stationary Wigner function distribution 
(Fig. 7.3) and verify the scaling given by the temperature dependent length 
L 




Figure 7.3: Tunneling Wigner distributions as a function of the temperature. The 
relevant parameters are: 7 = 0.25a;, T = O.OOlu;, n R = 0,0.75,1.5 respectively 
represented by dots, circles and asterisks. Full lines are the analytical results. See 
Fig: 7.1 for the other parameters and the scales. 



7.1.3 Current 

The current through the SDQS towards the right lead is given by the general 
formula: 

00 

I R {t) = Y, n Pn{t) (7.25) 

n=0 

In the tunneling limit the probabilities P n (t) are be expressed in terms of 
Pii\t) since: 

P n (t) = TYoscfpSWMO) +P { n ] (*KhK)] = P$(t) +P ( n\t) (7.26) 

Apart from the renormalized electronic rates the system has no sign of the 
oscillator degree of freedom and can be treated formally as a static quantum 
dot. We write in complete analogy with (3.54) the right (left) lead time 
dependent current: 

I R (t)=f RPll (t) 
h(t)=f L Po (t) 

116 



7.1. RENORMALIZED RESONANT TUNNELING 



0.18 
0.16 
0.14 
0.12 

■ « 

- 0.08 
0.06 
0.04 
0.02 




0.05 



Current 



0.1 0.15 

y / co 



Numerical 
Tunneling 



0.2 



0.25 



Figure 7.4: Current as a function of the damping for the SDQS. The asymptotic 
tunneling limit is indicated. The parameters are: Tl = Tr = O.Olu;, 7 = 0.25w, 
d = 0.5xq, A = 2xq, T = 0. 



In the stationary limit they coincide: 
lim^oo I R (t) = f R pf{ 

rstat / 

~ \ 

We stress the difference between Tl(r) and Tl(r)' while the first is the average 
injection (ejection rate) and thus represents the current through the device, 
the second is an average over the mechanical distribution function of the 
position dependent rate operator. Despite being an average rate it does not 
take into account the charge state of the dot (and the Coulomb blockade) and 
thus cannot represent the current through the device. We show in figure 7.4 
the current calculated numerically and the asymptotic value of the tunneling 
regime given by Eq. (7.28). 



R 



r R r 



RL L 



r 



(7.28) 



7.1.4 Current-noise 

The calculation of the current noise and Fano factor in the tunneling regime 
is interesting since it gives us the opportunity to explain the concepts and 
methods we introduced in section 5.3 still keeping a computational level 
analytically affordable due to the much lesser number of degrees of freedom 
of the effective model. 
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We start with the MacDonald formula for the zero frequency current 
noise: 



S(0) = t hr^ ± [ n 2 P n {t) - ( nPn{t)) 2 } (7.29) 



n=0 



n=0 



where P n (t) = p^ + Pii an d the equation of motion for are: 



(n) 



■ (n) 

Poo 



— 1 lPoo 1 RPii 



(7.30) 



We identify the effective Liouvillean C for the relevant variables pu(t) in the 
model: 



-r L 



(7.31) 



The evaluation of the different terms of the current-noise can be carried out 
by introducing the generating functions Fu(t; z) = ^2 n Pii \t)z n and following 
the same reasoning we indicated in section 5.3.2. The result is formally the 
same 7 : we express the Fano factor in terms of the stationary probabilities 

P% 



,sl '' lt and the pseudoinverse of the the Liouvillean QC 1 Q. 



1 - 



Jstat 



Q C -lQ [ U >-R) (Pof 

o o J ^ ^ V o o ) lpfi at 



o v R \ „^„(o r 



(7.32) 



where || • ||i is the sum of the elements of the vector. The projectors V and 
Q, built from their definition (see sec. 5.3.2), have the matrix form: 



Poo Poo \ 
Pll Pll / 



<r> — I roo roo I n — I ^oo (7 W\ 

'" I - i ? ~ " I -!..: pStat I • l'- JJ j 

For the explicit calculation of Eq. (7.32) we define the auxiliary vector E: 



Pii 
'Pii 



Soo 
En 



that satisfies the equation 



,/0 I /A f/'on : 

o oi 1 p s A at 



QC~ l Q 



Pll 



r R (p s A at ) 2 



i 

-i 



(7.34) 



(7.35) 



7 The perfect matching can be better appreciated in a more general super-operator 
formulation in which the current noise reads in both cases ({0\Ir — 2I R Q£~ 1 QI R \0)) . See 
for example [33] for details. 
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Figure 7.5: Current-noise as a function of the damping for the SDQS. The 
asymptotic tunneling limit is indicated. The parameters are: Tl = T R = O.OIuj, 
7 = 0.25^, d = 0.5x , A = 2x , T = 0. 



Together with the condition £ o + Sn = 8 we can calculate £ and by 
substitution into (7.32) find the Fano factor F: 



1 2T R E U _ 1 2T R (T L + T R ) T R Tl Tj + r| 

^ stat TlTr (f L + r R f (f L + f R y 



7.2 Shuttling: a classical transport regime 

The simplified model for the shuttling dynamics is based on the observation - 
extracted from the full description- that the system exhibits in this operating 
regime extremely low Fano factors (F ~ 10~ 2 ). In the simplified model we 
assume that there is no noise at all in the system. Its state is represented by 
a point that moves on a trajectory in the 3-dimensional device phase-space 
of position, momentum and charge of the oscillating dot. It is hard not to 
imagine the charge on the oscillating dot as a stochastic variable governed 
by tunnelling processes. Nevertheless in the shuttling regime the tunneling 
events are made effectively deterministic since they are very highly probable 
at specific times settled by the mechanical dynamics. On the other hand, 

8 This condition is readily proven since the last operator to the left in the definition of 
£ is Q and the sum on the elements of E would extract the component in the direction of 
the null vector of the Liouvillean. This component has been already projected out by Q. 
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due to the electrostatic force e£ acting on the charged dot that mirrors 
the electrical dynamics into a mechanical acceleration only "deterministic" 
tunneling events can produce a regular mechanical dynamics. We discuss in 
the next session the derivation of the equations of motion for this semiclassical 
model 9 and then compare the results with the full quantum description. 

7.2.1 Equation of motion for the relevant variables 

The dynamics of the SDQS is described by the set of coupled Klein-Kramers 
equations (4.32). In order to implement the zero noise assumption we first 
set T = and simplify the equations further by neglecting all the terms of 
the H expansion since we assume the classical action of the oscillator to be 
much larger than the Planck constant 10 . We obtain: 



dr 



x 9 
X dP 



p 



d_ 

dX 



uodP 



W, 



oo 



LO CO 



d\ d 



d 



X — — I — - P— J- 1 9 p 
A J OP dX udP 
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L -2X 



11 



w, 



cl 
00 



where we have introduced the adimensional variables: 



T = CUt, X = 



1 

A' 



P = 



V 



mwX 



(7.37) 



(7.38) 



The superscript "cl" indicates that we are dealing with the classical limit 
of the Wigner function justified by the complete elimination of the quan- 
tum "diffusive" terms from the Fokker-Planck equations and thus indefinitely 
sharp distribution functions are now allowed. We can at this point introduce 
a separation Ansatz: 



W/ C ' (X, P, r) = p o(r)5(X - X*(r))6(P - P cl (r)) 
Wf 1 (X,P,T)=p n (r)5(X - X d (r))5(P - P cl (r)) 

9 The derivation is not rigorous and the equations of motion (7.43) should be probably 
viewed as an "educated guess" to describe the shuttling device. Nevertheless the effort for 
a mathematical derivation from the Klein-Kramers equations (4.32) has been motivated 
by the very good agreement between the extremely simple model and the full description. 

10 This condition is fulfilled in the shuttling regime since the area enclosed by the ringlikc 
structure of the total Wigner distribution in the oscillator phase space is much larger than/l. 
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where the trace over the system phase-space sets P00 + P11 = 1- The variables 
X cl and P cl represent the position and momentum of the (center of mass) 
of the oscillating dot; pn(oo) is the probability for the quantum dot to be 
charged (empty). 

By inserting the Ansatz (7.39) into equation (7.37) summing over the 
charge index % and matching the coefficients of the different distribution func- 
tions on the left and right hand side we obtain: 



P- = _*- + ' Ip- < 7 - 40 > 

A uo 

In order to close the system of differential equations we need the equations 
of motion for p u 11 provided by the integration over X and P of (7.37) with 
the separation Ansatz: 

Poo = e poo H e p u 

r U r U (7.41) 

1 L -2X rA 1 R 2X cl 

Pn = — e poo e p u 

uo uo 

Unfortunately the system of equations (7.40) and (7.41) is not equivalent 
to (7.37): integration and summation have introduced spurious solutions. 
If we substitute the Ansatz (7.39) into the equations (7.37) and use the 
equations of motion (7.40) and (7.41) we obtain the condition: 

PooPu = (7.42) 

The only differentiable solution for this equation is poo = or pu = for 
all times and is not satisfying the equation of motion (7.41). Strictly speaking 
there are no solutions of (7.37) in the form suggested in the Ansatz 12 . 

Nevertheless we can imagine that the electrical switching time between 
the only allowed charged {pu = 1; Poo = 0) and empty (poo = 1; Pu = 0) 
states is much shorter than the shortest mechanical time (the oscillator period 
T = 2n/uo ). A solution of the system of equations (7.40) and (7.41) with 
this time scale separation could "almost everywhere" satisfy the condition 
(7.42) and, inserted into (7.39) represent a solution for (7.37). 

We rewrite the set of equations (7.40) and (7.41) as: 



11 The dynamics of poo and Pu are obviously coupled. 

12 First order partial differential equations can be solved by the characteristics method 
that calculate the flow along trajectories: this is not true for system of partial differential 
equations. More physically: the switching process in itself is introducing noise. 
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X = P 

P = -X + d*Q - 7 *P (7.43) 
Q = r L e~ 2X (l-Q)-T* R e 2X Q 

where we have dropped for simplicity the "cl" superscript from the me- 
chanical variables, we have renamed pn = Q and used the trace condition 
Poo = 1 — p\\. We have also defined the rescaled parameters: 

d* = d/\, t* = 7/w, Tl R = T LtR /uj (7.44) 
7.2.2 Stable limit cycles 

The set of coupled non-linear differential equations (7.43) represents the 
starting point of the analysis of the simplified model for the shuttling regime. 
We calculate the numerical solution for different values of the parameters and 
different initial condition. For the parameter values that correspond in the 
full description to the shuttling regime, the system has a limit cycle solution 
with the desirable time scale separation we discussed in the previous section. 
We report in figure 7.6 the typical appearance of the limit cycle. In the 
first graph (a) we show the charge Q(r) as a function of time. The charge 
value is jumping periodically from to 1 and back with a period equal to 
the mechanical period. The transition itself is almost instantaneous. The 
limit cycle is a trajectory in the 3D phase-space of the device (X,P,Q). In 
the graphs (b), (c) and (d) of figure 7.6 three different projections of the 
trajectory are reported. The X, P projection shows the characteristic cir- 
cular trajectory of harmonic oscillations. In the X, Q (P, Q) projection the 
position(momentum)-charge correlation is visible. From the combination of 
the two projection we can argue that the trajectory in the X, Q graph is 
drawn clockwise during the cycle. The oscillating dot gets charged on the 
left, it is then carrying the charge towards the right and finally it unloads 
the electron to the right lead before returning empty towards the left. The 
amplitude of the oscillation is several times the tunneling length. All the 
qualitative features of the shuttling regime are present. 

The full description of the SDQS in the shuttling regime has a phase space 
visualization in terms of a ring shaped stationary total Wigner distribution 
function. We can interpret this fuzzy ring as the probability distribution 
obtained from many different noisy realizations of (quasi) limit cycles. The 
stationary solution for the Wigner distribution is the result of a diffusive 
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dynamics on an effective "Mexican hat" potential 13 that involves both am- 
plitude and phase of the oscillations. In the noise-free semiclassical approxi- 
mation we turn off the diffusive processes and the point-like state describes 
in the shuttling regime a single trajectory with a definite constant amplitude 
and periodic phase. We expect this trajectory to be the average of the noisy 
trajectories represented by the Wigner distribution. In the third column 
of figure 7.7 the total Wigner function corresponding to different parame- 
ter realization of the shuttling regime is presented. The white circle is the 
semiclassical trajectory. In the first two columns the asymmetric sharing of 
the ring between the charged and empty states is also compared with the 
corresponding Q — 1 and Q = portions of the semiclassical trajectory. 

In the semiclassical description we also have direct access to the current 
as a function of the time. For example the right lead current reads: 

Ir(t) = Q(r)T R e 2X ^ (7.45) 

In figure 7.8 the right current is presented as a function of time for a few 
mechanical oscillation periods. The current is given by a series of spikes 
that well represent the single electron released to the right lead after being 
shuttled by the oscillating dot. 

The average current is a coarse grained measurement of the time depen- 
dent current. Since the current is periodic, if the measurement time is long 
enough we obtain a stationary current: 

P t&t = ^J T ^dr'I R {r') (7.46) 

Since the current is a periodic function of time the average over one period 
is equal to the average over an infinite time. The numerical integration of 
the function plotted in figure 7.8 gives 

jstat = 0.15916 (7.47) 

which is with impressive accuracy the "magic value" of 1/27T (i.e. one electron 
transferred per oscillation period). 

Also for the TDQS we tried a semiclassical analysis of the shuttling 
regime. For the derivation of the semiclassical equations we start again from 
the Hamiltonian (2.9) but we consider x and p as classical variables. The 
GME equation for the reduced density matrix a is derived by keeping x as 
a parameter (p does not enter explicitly the formulation of the GME). Since 
the oscillator is treated classically also the elements of the density matrix 

13 See the next section for details on the derivation of this effective potential. 
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Figure 7.6: Different representations of the limit cycle solution of the system of 
differential equations (7.43) that describes the shuttling regime, a) Charge on the 
dot as a function of time. The charging and discharging times are the shortest 
time scales, b) Circular trajectory in the mechanical phase space: the motion of 
the dot is harmonic, c) Projection of the limit cycle in the charge position plane: 
the trajectory shows charge-position correlation, d) Projection of the limit cycle in 
the charge momentum plane: the trajectory shows charge-momentum correlation. 
X is the coordinate in units of the tunneling length X, while P is the momentum 
in units ofmcoX. 
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Figure 7.7: Correspondence between the Wigner function representation and the 
simplified trajectory limit for the shuttling regime. The white ring is the (X, P ) 
projection of the limit cycle. The Q = 1 and Q = portions of the trajectory 
are visible in the charged and empty dot graphs respectively. The parameters are 
7 = 0.02a;, d = 0.5xo, T = 0.05w in the first two rows and V = 0.01 in the last 
one. The tunneling length is X = xo in the first row and A = 2xo in the second 
and the third. 
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Figure 7.8: Time resolved shuttling current calculated on the limit cycle solution 
of the system of differential equations (7.43). The regular spikes well represent the 
single electron shuttled per cycle in this regime. 



aij are now real numbers and the trace sum rule can be used to reduce the 
effective number of equation for the electrical dynamics 14 from 10 to 9. The 
GME in matrix notation reads: 

cr = - l -[H sys ,a] + E[a] (7.48) 

where ^ 

1 .2^2 , P , I j. /_\ AV, 



H sys = -mu 2 x 2 + + ( t L (x) -^x t R (x) I (7.49) 

t R (x) 



2 2m 
and the driving Liouvillean H has the form 



AV 
2 



(1 - o~ LL - o cc - a RR -0"lk/2\ 
-cr CR /2\ (7.50) 

-a RL /2 -a RC /2 -a RR J 

where we have taken = T R = T. As equation of motion for the mechanical 
variables x and p we take the Hamilton equations derived from the effective 
Hamiltonian 

(H Bys ) = Tr cl [(7# sys ] (7.51) 



14 Also in this case the elements <7oi an d °"i0 with i — L,C,R vanish identically reducing 
the effective number of equations from 16 to 10. 
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where Tr el is the trace on the electronic states, combined with a phenomeno- 
logical friction: 

± _ d(H sys ) = p_ 
dp m 

d(H sys ) 2 AV f 7 r 9 \ 

p = IV = —mw x + - — occ - IP v ' - oz ) 

ox 2xq 

+ y (Pcl + cr LC ) - y (<J C r + cr RC ) 

The phase space has now 11 dimensions (9 electrical and 2 mechanical) 
but the principle is the same as in the SDQS. Also in this case the solution of 
the system of differential equations (7.48) and (7.52) exhibits for parameters 
that correspond to the shuttling regime a stable limit cycle solution. In figure 
7.9 we compare the Wigner distribution function for the central dot with the 
projection on the (x,p) plane of the limit cycle trajectory of the semiclassical 
approximation. The probability to find the electron on the central dot ace 
is not in general oscillating between and 1. The red dots in the figure 
indicate the position of the maximum occupation. In the TDQS the half 
moon characteristic of the shuttling regime rotates in the phase space [10], 
close to the first two electromechanical resonances, as a function of the device 
bias. The correspondence between the darker areas of the Wigner function 
and the red dots in figure 7.9 shows that also this property is captured by 
the semiclassical approximation. 
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Figure 7.9: Comparison between the semiclassical approximation and the full de- 
scription for TDQS. The figure shows the central dot Wigner function and the 
respective limit cycle. The dark dot on the trajectory is in correspondence with the 
maximum of the semiclassical occupation ace- Also the phase space rotation is 
captured by the simplified model. The matching is best in the shuttling regime. See 
figure 5.4 for the details on the parameters. 
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Figure 7.10: Schematic representation of the time evolution of the current in the 
dichotomous process between current modes. The scheme refers in particular to 
the SDQS coexistence regime where the currents involved are the shuttling (I s h) 
and the tunneling (Itun) currents and the names of switching rates T[ n and T ou t 
indicate the behaviour of the dot oscillation amplitude (the relevant variable in the 
model) during the transition. 



7.3 Coexistence: a dichotomous process 

The slowest dynamics in the coexistence regime is represented by the switch- 
ing process between the shuttling and the tunneling regime. The amplitude 
of the dot oscillations is the relevant variable that is recording this slow dy- 
namics. We analyze this particular operating regime of the SDQS in three 
steps. We first assume the slow switching mode between two different cur- 
rent channels and explore the consequences of this hypothesis in terms of 
current and current noise. We then derive the effective bistable potential for 
the amplitude associated with the slowest dynamics of the shuttle device in 
the coexistence regime. Finally we apply Kramers' theory for escape rates 
to this effective potential and calculate the switching rates between the two 
amplitude equilibrium states corresponding to the local minima of the po- 
tential. We conclude the section comparing the (semi) analytical results of 
the simplified model with the numerical calculations corresponding to the 
full description. 

7.3.1 Dichotomous process between current modes 

Let us consider a bistable system with two different modes + and — and two 
different currents I + and I_, respectively, associated with the modes. The 
system can switch between mode + and — randomly, but with definite rates: 

namely a for the process H > — and j3 for the opposite ► +. We collect 

this information in the master equation: 
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d_ fP + \ f-a /3\ (P, 
dt \P- \ a -P) VP_ 



LP (7.53) 



where P + and P_ are the probabilities to find the system in the mode "+" 
and "— " respectively. 

The matrix notation suggests to read the equation (7.53) at a more general 
level and embed its classical rate dynamics into a quantum formalism. 

Instead of the modes + and — we consider the state vectors |+) and 
|— ): they define the Hilbert space of the quantum system. Density operators 
defined on this Hilbert space represent the state of the system. Since we 
are only interested in classical states we assume that the density matrix a 
is always diagonal in the |+), |— ) basis. In terms of the probabilities P the 
density matrix reads obviously: 

_ fa++ a+\ (P + 



, - \ n PI ( 7 - 54 ) 

Since each of the two states is associated to a well defined current, also 
the current operator is diagonal in the +, — basis: 



I + 
/_ 



(7.55) 



The average current is the trace of the current operator multiplied by the 
density matrix, namely: 

(/(*)> = Tr[a(t)I] = P + (t)I+ + P-(t)I- (7.56) 

and it is in general time dependent since the occupation probabilities evolve 
in time according to the equation (7.53). 

Equation (7.53) is the equation of motion for the density matrix a. In 
general the Liouville space 15 for a system with a Hilbert space of dimension 
N has dimension iV 2 . In this particular case the system is simplified: the 
space of the diagonal linear operators has again dimension 2. The equation 
(7.53) is the representation of the equation: 



in the basis: 



a = C[a\ (7.57) 



I,', u) - = (S f) 



with the identification: 



15 



It is the space of the linear operators on the original Hilbert space. 
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<-p=(£), r^L=(- J) ,,59, 

Until now we have used the Schrodinger picture in which the states evolve. 
For the definition of the current noise we have to introduce the analogue 
for open irreversible systems of the Heisenberg picture defined for closed 
reversible system 16 . The equation (7.57) has the formal solution for a time 
independent Liouvillean: 

a(t) = e tc a(0) (7.60) 

The initial condition <r(0) is evolved to the time t by the super-operator 
e tc . The Heisenberg picture for the operators is designed to preserve the 
average: 



(I(t)) = Tr[I s a(t)] = Tr[Ise tC a(0)} = Tr[e tC ' l s a(0)] = Tr[I H (t)a(0)} (7.61) 

where the (current) operator appears first in the Schrodinger Is, and then 
in the Heisenberg picture. & is the adjoint of the superoperator C. It 

is useful to consider the representation of the above equations in the basis 
(7.58): 



(I(t)) = (I S , P(0) = (I, e tL P(0)) = (e* Lt I, P(0)) = (I H (t), P(0)) (7.62) 

where Is = [I+, I-] T is the vector representation of the current operator in 
Schrodinger picture. It is also not difficult to realize that the trace of the 
product of two operators is equivalent to the scalar product in the vector rep- 
resentation. The current-noise is the spectral density of the average current 
fluctuations: 



/+oo 
dr[(I(t + r) J(f )) - (I(t + r)> (/(*)>] e wr (7.63) 
-oo 

In terms of the current operators in the Heisenberg picture, the noise reads: 

/+oo 
dr{Tr[I H (t + r)I H (t)a(0)} 
oo (7.64) 

-Tr[I H (t + r)a(0)}Tr[I H (t)a(0)]}e^ T 



16 



For a rigorous treatment of the argument sec for example H.-P. Breuer et al. [46]. 
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The system tends for large times to the stationary state: 



jstat 



jstat 



g 

a + (3 
a 

a + 13 



(7.65) 



as can be verified by substitution in (7.53). The corresponding stationary 
current reads: 

jstat = Um y( t y = hP + I-a (7 g6) 



a + (3 



We are interested in the zero frequency noise spectrum uo — > 0. We split 
the integral into two parts: over negative and positive r respectively. Using 
the symmetry under time reversal r — > — r of the integrand we obtain: 



5(0) 



lim lim 23? 

U>^0 t^OO 



rfr fTr[/^(t + r)/ H (t)<r(0)] - / stat2 J e t{uJ+t£)r 



(7.67) 

where we have also added to the frequency the small imaginary convergence 
factor is. We now use the quantum regression theorem to rewrite the current- 
current correlator and perform the limit in t: 



S(0) = lim 23? 

LU—*0 



dr (Tr[/e r£ /(7 stat ] - J stat2 ) e i(w+ie)T 



We then integrate in r and obtain: 



(7.68) 



S(0) = lim 23? 



Tr[/(£ + jw)-70 - 



(7.69) 



In the limit u> — * the resolvent (£ + iuo)' 1 is singular since £ has a non- 
trivial null space. To handle this singularity we introduce the super-operator 
projectors: 

p[m] = (j stat Tr[«], Q = l-V (7.70) 

We already encountered these projectors in section 5.3.2 where we discussed 
their properties (5.34). In terms of these projectors the resolvent can be 
rewritten: 



(C + iu)- 1 = Q(C + iuu)- l Q + — 

too 



(7.71) 
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The term Q(C + iuj)~ l Q is well-behaved in the limit to — > since the Q 
projectors are confining the inverse of the Liouvillean on the subspace or- 
thogonal to the null space. The term ?j is still singular. However, using the 
definition (7.70) we can calculate: 



Tr 



V 

I—Ia st&t 


= Tr 


rstat 
j ^.stat 


iu 




IUJ 



rstat 2 



IUJ 



(7.72) 



and notice that the two diverging contribution in the current noise exactly 
cancel. The zero frequency current noise reads: 



5(0) = -2Tr(/Q£- 1 Q/a stat ) 



(7.73) 



For the evaluation of the current noise 5(0) we introduce the auxiliary 
quantity £ 

E = QC- l QIa st&t (7.74) 

so that 

5(0) = -2Tr(/S) = -2(1, S) (7.75) 

where we have used in the last equality the vector notation characteristic of 
the example at hand. The quantity £ satisfies the equation: 

££ = QIa stat = (I - / stat ) (7 stat 



that in vector notation reads 



-a (3 
a —f3 



(j± ~ I- 

(3 



a 



a 



-0 



(7.76) 



(7.77) 



The components of the vector S also satisfy the independent relation S + + 
£_ = 0. We solve then for S and we are able to give an analytical expression 
for the noise of the dichotomous process: 



5(0) 



af3 



:A/ 2 



(7.78) 



where AI = — Summarizing, we give the expression for the stationary 
current and Fano factor for a dichotomous process between two current modes 
with currents I + and 7_ and switching rates a and j3: 



rstat 



I+/3 + I-a 
a + 13 
5(0) a(3 



(A 



(7.79) 



J stat (a + (3) 2 I +13 + I _a 
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Figure 7.11: Schematic representation of the effective potentials for the three 
operating regimes. 



The framework of the coexistence simplified model is given by these re- 
sults. The task is now to recognize in the dynamics of the shuttle device the 
two modes and above all calculate the switching rates. The answer are in 
the Kramers' escape rates for a bistable effective potential. 



7.3.2 Effective potential 

The stationary total Wigner function for the SDQS evolves as a function of 
the mechanical damping from a fuzzy dot close to the origin of the phase space 
(tunneling regime) to a growing ring-shape in the small damping (shuttling) 
regime. At intermediate damping rates the shuttling ring and the tunneling 
fuzzy dot coexist. Those properties of the Wigner function can be understood 
in terms of an effective stationary potential in the phase space generated 
by the non-linear dynamics of the shuttle device. We show in figure 7.11 
the three qualitatively different shapes of the potential guessed from the 
observation of the stationary Wigner functions associated with the three 
operating regimes. Fedorets et al. [9] started the analytical work for the 
understanding of the tunneling-shuttling transition 17 in terms of an effective 
radial potential. Taking inspiration from their work we extend the analysis 
to the slowest dynamics in the device and use quantitatively the idea of the 
effective potential for the description of the coexistence regime. 



17 Very recently they also introduced the spin degree of freedom in the oscillating dot 
[47]. 
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Elimination of the fast dynamics 

The starting point are the Klein-Kramers equations for the SDQS that we 
rewrite symmetrized by shifting the coordinates origin to d/2: 
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[7.80) 



One by one we will now get rid of the variables that due to their fast dynamics 
are not relevant for the description of the coexistence regime. In equations 
(7.80) we describe the electrical state of the dot as empty or charged. We 
shift to the description in terms of state + and state — with the definition: 

W + = W 00 + W U , W- = Woo-W u (7.81) 
In terms of these new Wigner functions the Klein-Kramers equations read: 



d t W + = {H osc , W + } P + C 1 W + + C 2 W- 

d t W_ = {H osc , W_} P + (A + T + )W^ + (C 2 + Y^)W + 

where 

T + = T R e 2q/x + T L e' 2q/x 
T_ = T R e 2q/x - T L e- 2q/x 

(1 \ F i / h \ 2n 

%+ 2 ^f^2d U) ^ (7 ' 83) 
' n=i ' ^ ' 




The symbol {•, «}p represents the Poisson brackets and we have introduced 
for the partial derivatives the short notation d x = J^. It is useful for the 
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calculation of the effective potential to treat position and momentum on 
equal footing. For this reason we introduce the adimensional variables: 

The general structure of the equations (7.82) remains unchanged but the 
Liouvilleans and the T functions take the form: 



1 + = — e H e 

to to 

T _ = Ev x _ T _L e -2x 

UJ UJ 

x ' n = ~\ 



n=l 

d a + r -V (_1)n P°W 



}2n 

2n! V A y L 

n=l 



Since in the following we will assume the "classical" limit of a tunneling 
length much larger than the zero point uncertainty in the position, only the 
first term of the sums that appears in the definition of C\ and £ 2 will be 
kept. Following Fedorets et al. [9] we also introduce polar coordinates: 

X = A sin 0, P = A cos (7.86) 
The equations (7.82) are transformed into: 



d T W + = (-fy + d)W+ + C 2 W_ 

d T w„ = (-fy + £i - r + )H/_ + (£ 2 + r_)H/+ 



(7.87) 



It is interesting to note that the Poisson brackets for the harmonic oscillator 
reduce in the polar coordinates to — d$ since there is no amplitude dynamics 
in the phase space. The T functions now read: 



£ " (7.88) 

Y _ zJl_ e 2Asm<f> _ 1 L 2Asm<f> 
UJ UJ 

For the transformation of the Liouvilleans C\p we used the polar coordinate 
representation of the following differential operators: 
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<9 P = — — y + cos 0^ 

=— cos 4>dAA — -dcj) sin 

./i ./i 

02 2 . sin 2 sin 2 o2 

9 



sin0cos09 A0 + cos 0<9 A ^ gg ^ 

1 2 1 

=- cos 2 09iA - -d 2 M sin cos + — d\ sin 2 

— —(9^4 sin 2 — tj<9</> cos 4> sm 

dpP =— cos 2 0<9yiA 2 — d<f, sin cos 
= 1 — sin cos 0(9^ + cos 2 $A3a 

where d 2 A<t> = qj^- The two different forms (partial differential operator &a 
and/or to the extreme left or right) that for the moment seem redundant 
will be very useful in the derivation of the effective potential. 

We now start to make approximations on the system of equations (7.87) 
with the main idea of eliminating the fast irrelevant variables. We know that 
eventually the density operators an that describe the system will reach a 
stationary state. In absence of the harmonic oscillator the state + would be 
fixed by the trace sum rule and the state — would relax to zero on a time 
scale fixed by the tunneling rates. We assume that also in the presence of 
the mechanical degree of freedom the relaxation dynamics of the |— ) state 
is much faster than the one of the |+) state. We set d T W^ to zero in the 
second equation of (7.87) and formally solve the equation for W_. We insert 
the result in the first equation of (7.87) and obtain a closed equation for the 
(total) Wigner function W + : 



d T W + = [-fy + A + £ 2 (1 - GoA)- 1 ^/^ + T_)]W + 

" v ' (7.90) 

where Go = (<9</> + T + ) _1 and we also define the (super)operator C + . At this 
point the detailed information about the charge state of the oscillating dot is 
integrated out, but the equation (7.90) still takes the fast phase dynamics into 
account. Since we are interested only in the amplitude dynamics (the slowest 
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in the coexistence regime) we introduce the projector V4, that averages over 
the phase and its orthogonal complement Q4, — 1 — V4,: 

1 f 27T 

nM = ^ J d<t>. (7.91) 

Using these two operators we decompose the Wigner distribution function 
into: 



W+ = VW+ + QW+ = W+ + W+ (7.92) 

We apply the same projectors also to the equation (7.90) and we obtain the 
set of coupled partial (integro- 18 ) differential equations: 



d T W + = V^C + W + + p*c + w + 

d T W + = [-fy + Q £ + ] W + + Q 4 ,C + W + 



where we have used the properties: 

V+dt = 



Qthdd, = 84,0.4, = 84 



(7.94) 



which is readily demonstrated since all the distributions involved are periodic 
functions of 0. We now assume that the phase relaxation time that governs 
the dynamics of the phase-dependent distribution W + is much shorter than 
the amplitude relaxation time of the distribution W + 19 . We then set d T W+ = 
in the second equation of (7.93), solve it for W + and insert the result in 
the first equation. The result is a closed equation for the (quasi)probability 
distribution function W<: 



d T W + = V4,C + [l + (1 - g Q4,C + )- 1 g Q^+W + (7.95) 

where go = ((9^) _1 . We notice that the operator g Q is in the equations always 
following the projector Q4. This combination of operators ensures us as an 
output a periodic function in the variable 0. 

18 The operator C contains the inverse of differential operators. 

19 This assumption is based on the observation that the diffusion dynamics that governs 
the relaxation times is largely facilitated along the transverse direction where no potential 
barriers (see fig. 7.11) must be overcome. The potential barriers occur instead in the 
radial direction. 
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Small parameters expansion 

We believe that equation (7.95) contains the relevant dynamics of the SDQS 
in the coexistence regime, but, despite all the approximations already in- 
troduced the problem still looks untractable. Again following the already 
mentioned work by Fedorets et al. [9] we consider the perturbation expan- 
sion of (7.95) in the "small parameters": 



These three inequalities describe a limit often encountered in this thesis and 
they correspond respectively to the three physical assumptions: 

1. The external electrostatic force is a small perturbation of the harmonic 
oscillator restoring force in terms of the sensitivity to displacement of 
the tunneling rates. This ensures a quasi oscillator-independent treat- 
ment of the tunneling regime. 

2. The tunneling length is big compared to the zero point fluctuations. 
Since the oscillator dynamics for the shuttling (and then partially also 
for the coexistence) regime happens on the scale of the tunneling length, 
this condition ensures a quasi-classical behaviour of the harmonic os- 
cillator. 

3. The coupling of the oscillator to the thermal bath is weak and the 
oscillator dynamics is under-damped. 

We want to expand the equation (7.95) up to second order in the three 
small parameters (7.96). To proceed systematically we first notice that they 
appear at least to first order in the Liouvilleans C\ and £ 2 and thus also in 
the operator C + . We can then safely expand (7.95) up to second order in C + 
without missing any desired term: 



Using the definition of C + contained in (7.90) we now expand this last equa- 
tion up to second order in C\ and L<i- 




(7.96) 



(7.97) 




(7.98) 
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The last step is to consider within the equation (7.98) only the contributions 
up to second order in the expansion parameters (7.96). In the definition of 
the Liouvilleans in the adimensional form (7.85) we note that each of the 
small parameters (7.96) is associated with a differential operator dp. We 
then average out the phase variable from the equations. We thus expect the 
second order expansion in the small parameters (7.96) to be of second order 
in the differential operator 8a- More precisely the expansion takes the form: 



where V'(A) = j^V(A) and D(A) are given functions of A. Before calculating 
explicitly all the different contributions that compose the functions V and D 
we want to explore the consequences of the formulation of the Klein-Kramers 
equations (7.80) in the form (7.99). The stationary solution of the equation 
(7.99) reads: 



where Z is the normalization that ensures the integral of the phase-space 
distribution to be unity: dA'2nA'Wf at (A') = 1. The equation (7.99) is 
identical to the Fokker-Planck equation for a particle in the bidimensional 
rotationally invariant potential V (see figure 7.11) with stochastic forces de- 
scribed by the (position dependent) diffusion coefficient D. 

All the differential operators acting on W+ in (7.98) are at maximum of 
second order in 8a- In order to obtain Eq. (7.99) we substitute into (7.98) the 
definition of C\ and £2 and we "push" the differential operator 8a all the way 
towards the left or the right in the sequence of operators acting on W+. All 
the components in which we are left with only one 8a operator on the right 
define the effective potential V. All the others have the second 8a on the 
extreme right and define D(A). All contributions to the effective potential 
V and diffusion coefficient D can be grouped according to the power of the 
small parameters that they contain. 



8 T W + (A,r) = jd A A\V'(A) + D(A)d A ]W + (A,T) 



(7.99) 




(7.100) 



(7.101) 
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where the a functions read: 



«o = P<j> cos0G o r_ 



ai = -^cos0r_9 P (G o r_) 
«2 =7^cos0(%r_# o Q0dp(G o r_) 
a 3 =P cos0[G^r_ + ^^(Gor-) 



-sin09p(G o r_) 



(7.102) 



and the /5's can be written as: 



A =^cos 2 0[r + -r_G o r_ 




- sin 20Gor 



(7.103) 



In the expression above the variable disappears on the RHS when we ap- 
ply the projector V<f>. We have kept the somehow mixed notation with the 
differential operator dp to keep the notation lighter. In general the operators 
are acting on all what is on their right. Otherwise parentheses are limiting 
their range of action (see for example dp and in a\). As examples of the 
arguments that appear in the calculation of the a and j3 functions we give 
first the derivation of a "missing" contribution, then we explicitly derive 0:3 
and (3 3 . 

Some second order contributions in the parameters (7.96) simply do not 
appear in the final expansions 20 because the small parameter combinations 
which they represent is not present. The contribution in (^) instead is 
formally there but vanishes identically. From the expansion in terms of the 
operators C\ and £2 this contribution reads: 



We insert the polar coordinate expressions for the differential operators dpP 
(see Eq. 7.89) choosing the left operator form for the first and the right form 
for the second and we obtain: 



where we have used also the fact that the operator is applied to W + which 
does not depend on <f>. Remembering that = l—V^ and g is the indefinite 

20 It is the case of the contributions in (x) 2 an d 4\ (ir) 2 - 



'77 — 



V^dpPgoQ^dpP 



(7.104) 



'77 — 




(7.105) 
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integral in the variable we conclude the calculation: 
C 



-77 



^d A A 2 V^ cos 2 <f)g (cos 2 - - ) Ad A 



1 

2A 



(7.106) 



d A A 2 V^ cos 3 sin <pAd A = 



The functions a% and fa are both derived from the contribution of 
the small parameter expansion (7.98). The corresponding Liouvillean reads: 



C ld =V (j> dpG dpPG r.+ 
V lP d P Pg Q^d P G r_+ 
V^dpG T_g Q^d P P 



(7.107) 



We now use polar coordinates and take into account that the Liouvillean is 
applied to a function W + independent of the variable 0. We obtain: 



£ 7C z = —dj^AVtp cos 



G 2 T_ + G A cos 0<9 P (G o r-) + G Acos 2 <pG T_d A + 
A cos (t)goQ<t>dp(G T_) + Acos0# o Q</> cos 0G , o r_<9 A + 
GoT^goQcpA cos 2 <pd A 



The a 3 and (3 3 contributions can be split to obtain 21 : 



(7.108) 



C- d =- A d A A |^cos0[G' 2 r_ + AG d P (G T_) - ^ sin 4>d P {G Q V. 



A 
1 



d A A < Vs cos < 



+ 



G cos 2 0G o r- + -G r_ sin 20 - - sin 20G o r- 



(7.109) 

Since we have projected out the phase we are effectively working in a one- 
dimensional phase space given by the amplitude A. Equation (7.99) though is 
not as it is a Kramers equation for a single variable. This is related to the fact 
that also the distribution W+ is not the amplitude distribution function, but, 



21 In this passage we have also used the projector to define a scalar product 
TWWsW = if, 9) and the adjoint relation: (/, 6g) = (6 f f,g). 
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Figure 7.12: Example of the potential V (full) and the geometrically corrected po- 
tential V (dashed) for the SDQS in the coexistence regime. Notice for the corrected 
potential the logarithmic divergence in the origin. 



so to speak, a cut at fixed phase of a two dimensional rotationally invariant 
distribution. The difference is a geometrical factor A that it turns out is also 
simplifying the equation. We define the amplitude probability distribution 
W(A,t) = AW + (A,t) and insert this definition in equation (7.99). We 
obtain: 



which for an amplitude independent diffusion coefficient gives a corrected 
potential logarithmically divergent in the origin. The integration extremal 
is arbitrary and reflects the arbitrary constant in the definition of the po- 
tential. The equation (7.110) is the one- dimensional Kramers equation that 
constitutes the starting point for the calculation of the switching rates that 
characterize the coexistence regime. 

Evaluation of the and (3i functions 

We evaluate the «j and (3i functions numerically. Fortunately there is a 
rather natural discretization of the operators that compose the a and /3. 



d T W(A, r) = d A A[V'{A) + D{A)d A )jW(A, r) 
= d A [V(A)+D(A)d A ]W(A,r) 



(7.110) 



where we have defined the geometrically corrected potential 
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This discretization is given by the operator cos that opens the operator 
string of all the components and that can be considered as a sort of Fourier 
transform of the function that follows. We specify better the concept and 
the method by reducing explicitly to the "numerical form" the contributions 
a>o and a>2~- 



a =V<t> cos 0G o r__ ^ ^ 

«2 =V (f) cos(f)G T_g Q (f) dp(G T_) 

where we remind Go = {d<t> + F+)" 1 . We define the function G(A,<f>) = 
GoT-(A,(j)). If G n (A) is the discrete Fourier transform of G (periodic func- 
tion of 0) 

1 f 2lT 

G n(A) = —J o G(A,<p)e m * (7.113) 

with nfZ, then we can easily recognize the Fourier transform in the defini- 
tion of a and write: 

a (A) = ±(G 1 + G- 1 ) (7.114) 

We still need to identify the structure of the functions G n (A). It is useful for 
this purpose to calculate the Fourier transform of the functions T + and T_ 
and define the vectors f + and f _ whose components are: 

(F+) n = -[J n (-2iA) + J n (2iA)] 

% (7.115) 

(t.) n = -[J n (-2iA)-J n (2iA)] 
to 

where we have assumed for simplicity that = Tr = F and the Bessel 
function J n (z) is defined by: 

J n ( z ) = J- / e "^+ l2sin< ^0 (7.116) 

At the operator level we can write: 



o 



F_(A, 0) = [d* + F + (A, <f>)]G{A, 0) (7.117) 
We take the Fourier transform in on both sides and we obtain: 

oo _ 

(f _) n = y~] [in5 mn + (f +) nm ]G m (7.118) 

m =-oo 

lvl nm 
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where we have defined the matrix (f + ) nm = (r + )„_ m and the matrix M. Fi- 
nally we also define the Fourier transform of the cosine C = [. . . , \ , 0, \ , . . ] T 
with the 1/2 in position ±1 and we write the ao function in the compact 
numerical form: 



a = (C,M- 1 f_) (7.119) 

where the symbol (•, •) indicates the scalar product between the two 
vectors. All the vectors must be truncated in the numerical evaluation, but 
fortunately the Bessel functions J n decay fast with n and typically we verified 
numerical convergence taking only 40 terms around the zero position of the 
discrete Fourier transform "momentum space" . 

If the derivation of the numerical form for ao is revealing the spirit of the 
calculation, the a 2 case contains instead some of the typical "tricks". We 
start from the definition (7.112) and express the differential operator dp in 
polar coordinates. We obtain: 



a 2 = 7> cos0G o r-$oQ0 ^-^fy + cos^ G F_ (7.120) 
We absorb the partial derivative 8$ with the identification: 



d <t> = G 1 -T + (7.121) 

Another important element is the combination goQ^. It is important to study 
the Fourier transform of such an operator. Since g = (c^) -1 we obtain: 

i 

(go)mn = S nm (7.122) 

n 

The singularity for n = is cured by the projector that removes the 
th component in the Fourier transform of the vector on which it is acting. 
We call for this reason the combination goQ^ the pseudoinverse 22 of the 
differential operator d<f> and we use the symbol (d^ 1 )ps- 

The last problem that must be faced is the differential operator 8a that 
we always encounter in the expression OaG. Using the definition of G and 
the results obtained in the evaluation of a we can write: 

d A G = M^Oa? _ - M~ 1 d A (M)M- 1 f _ (7.123) 

22 From the arguments above is clear that goQ<f> = Q<f>goQ<p- 
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For the matrix M and the vector f _ we have the explicit formulation in terms 
of the Bessel functions. Using the recursive relation between the functions J 
and their derivatives 



we obtain: 

{d A M) mr 



d z J n (z) = -[J n -i(z) - J n +i(z)} 



iF 



to 



[J n _ m+1 (-2L4) - J n _ m _!(-2a) 
(+2iA) + J n _ m _!(+2a)] 



— [J n+1 (-2iA) - J n _!(-2iA) 
to 

+J n+1 (+2iA) - J n _!(+2L4)] = (<9 A f _) n0 
The "numerical form" of the function a 2 reads: 



(7.124) 



(7.125) 



a 2 = ~j(c, M _1 f_ (d<p)ps (5 f_ - Sf+G 
- A(7M _1 (9Af _ - d A (M)G) 

where we have also introduced the sin and cos Fourier matrices: 



(7.126) 
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(7.127) 



The other a and j3 functions are evaluated in a similar way and the calcula- 
tions involved do not present any further complication. 



7.3.3 Switching Rates 

We have proven in section 7.3.2 that the SDQS dynamics in the coexistence 
regime can be described by a Kramers equation for the amplitude probability 
distribution W with a given potential V(A) and diffusion "constant" D(A). 
We then dedicated the last section to describe an affordable and reliable 
numerical evaluation of the a and f3 functions that appear in the definition of 
the potential and diffusion constant. We are now able to identify completely 
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Figure 7.13: Example of the stationary distribution Wf 1 ^ (full) and the amplitude 
distribution yy stat (dashed) for the SDQS in the coexistence regime. The tunneling 
and shuttling states are in both cases well separated. 



the SDQS operating in the coexistence regime as a particular realization of 
the model for a dichotomous process between current modes presented in 
section 7.3.1. The shuttling and tunneling regimes with their characteristic 
currents are the two modes. Within the framework of the Kramers escape 
time theory we can now calculate the switching rates between these two 
modes. 

The effective potential V that we obtained has for parameters that corre- 
spond to the coexistence regime, a typical double well bistable shape (Fig. 7.12) 
We assume for a while the diffusion constant to be independent of the ampli- 
tude A. In this approximation the stationary solution of the equation (7.110) 
reads: 



where Z is the normalization Z = J °° W sta,t (A)dA. The probability distribu- 
tion is concentrated around the minima of the potential and has a minimum 
in correspondence of the potential barrier. If this potential barrier is high 
enough (i.e. V max — V m i n ^> D) we clearly identify two distinct states with 
definite average amplitude: the lower amplitude state corresponding to the 
tunneling regime and the higher to the shuttling (see fig. 7.13). 

The Kramers equation (7.110) describes the coexistence regime of a SDQS 




(7.128) 
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by mapping it into a classical model for a particle moving in a bistable po- 
tential V with random forces described by the diffusion constant D. It is 
this random force that allows the particle to overcome the barrier and jump 
between the two minima with definite rates. These rates are in the model 
for the SDQS the switching rates between the tunneling and shuttling mode 
Tin and r out that we introduced in section 6.3. The bistable potential model 
and the problem of calculating the average escape time (i.e. the average time 
necessary to leave one potential well) has been object of intense study be- 
cause, despite its simplicity, it finds application, in different local variations, 
in many branches of science 23 . The first result on this problem, published in 
1899 is the Arrhenius law for the chemical reaction rate k (which Arrhenius 
attributed to van't Hoff ) 

K = fle- Ec/kBT (7.129) 

where Eq is the activation energy, ks the Boltzmann constant, T is the 
absolute temperature and Q the attempt frequency or stearic factor. We 
base our derivation of the escape rates Tin and r ou t on the general treatment 
given in the book "The Fokker-Planck Equation Methods of Solution and 
Applications" by H. Risken [50]. 

Mean First Passage Time (MFPT) for a random variable 

Given a stochastic variable (in our case the amplitude A(t)) we can 
ask the question when this variable first leaves a certain domain (e.g. one of 
the potential wells). This time is called first-passage time. For definiteness 
we choose the initial condition for the amplitude A to be in the tunneling 
(lower amplitude) well and we ask when the amplitude is leaving for the first 
time the tunneling well to pass into the shuttling well. Since A is a random 
variable, also the first-passage time is a random variable since it depends on 
the specific walk of the amplitude. We want to calculate the distribution 
function for the first-passage times and in particular its first moment, the 
mean-first passage time. This time would constitute in our specific case 
the inverse of the tunneling to shuttling switching rate r out . We set the 
upper border A out between the saddle point amplitude A s and the shuttling 
amplitude A shut . We also set an arbitrary lower reflecting border in A min . 
24 (see Fig. 7.14). 

We define the probability density P(x, t\A tnn , 0) for the stochastic variable 
A(t) starting at t — with A(0) = A tun to be in x at time t without reaching 

23 For a review on this problem- also called the exit problem or the escape time problem- 
see for example [48] or [49]. 

24 The exact position of the borders is not relevant as far as they are far from the 
minimum or the maximum of the potential. 
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Figure 7.14: Bistable effective potential for the SDQS coexistence regime. The 
important amplitudes for the calculation of the rates are indicated. The reflecting 
(full) and absorbing (dashed) borders for the calculation of the escape rates. 



the upper border A ont . It is possible to demonstrate [50] that the probability 
P must obey the Kramers equation (7.110) if A min < x < A out , is identically 
zero for x > A out since the walks that touch the upper border are not taken 
into account and satisfy reflecting boundary conditions at the lower border, 
namely: 



— = C K P; P(x, 0| Aun, 0) = 6(x- Aun) for A min < x < A out 
P(x, t\A tun , 0) = for x = A out 
d x e^P(x,t\A tVLll ,0)\ x =A mta = 

(7.130) 

where Ck is the Kramers' Liouvillean of equation (7.110). The probability 
W(-^tun) t) of realizations which have started at A tnn and which have not yet 
reached the upper boundary up to the time t is given by: 

rAout 

W(A tun ,t)= / P( Xj t\A tun ,0)dx (7.131) 

The probability —dW of those realizations which reach the upper boundary 
in the time interval (t, t + dt) thus reads 

rAout 

-dW(A tw ,t) = - P(x,t\A tun ,0)dxdt (7.132) 

-^min 
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The distribution w(T) function for the first passage time T is therefore given 
by 



w(A tun ,T) = _ dW ^ T ) = _ [ A ° Ut P( x ,T\A tUQ ,0)dx (7.133) 



tun ) 

(7.136) 



dT 

" J 'mm 

The moments of the first-passage time distribution are 

poo pAout 

T n (A tun )= / T n w(A tun ,T)dT = / p n (x,A tun )dx (7.134) 

Jo J A min 

where p n (x, A tun ) is defined by 

PCX) 

p n (x,A tun ) = - / T n P(x,T\A tmi ,0)dT (7.135) 
Jo 

We obviously have 

poo 

Po(x, A tun ) = - / P(x, T\A tun , 0)dT = P(x, 0|Aun, 0) — S(x — A 
Jo 

Performing a partial integration gives the relation 

poo 

p n (x,A tun ) = n / T n - 1 P(x,T\A tun ,0)dT (7.137) 

Applying the Kramers' Liouvillean Ck to equation (7.137) we obtain the 
recursive relation: 

C K [p n (x,A tun )] = -np n -i(x, A tun ) (7.138) 

and in particular 

C K \pi(x, A tun )] = -5{x - Ann) (7.139) 

In other words the function pi(x, A tnn ) that enters the definition of the av- 
erage first-passage time is the Green's function of the Kramers' Liouvillean 
Lk- The Kramers' Liouvillean can be written in the form: 

V(A) . V(A) 

C K = Dd A e- — d A e + — (7.140) 

as can be easily proven by acting with the last partial derivative on the 
exponential. It is also possible 25 to give an analytic expression for the Green's 
function pi(x,A tun ) of the Kramers' Liouvillean Ck- 



25 



Basically it is a double integration of equation (7.139). 
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Pi(A,A tun ) 



D 



V(A) f^ out V(B)_ f 

d I e d 

J A IJA 



5(C - A tmi )dC 



dB (7.141) 



The mean-escape time from the tunneling well reads: 

i Mout r V(A) Mout V(B) \ r B ~\ \ 

T 1 (A tw ) = - \ e ~~ e— / 6{C-A tun )dC dB \ dA 

(7.142) 

The escape time from the shuttling well can be calculated exactly in the 
same way using A shut as starting point for the random process and A in and 
Anax as boundaries (see Fig. 7.14). The calculation of the integrals in (7.142) 
can be simplified considering that the integrand has a sharp maximum in the 
region (A, B) m (Ami, A) due to the behaviour of the effective potential V 
around those two points. The escape time from the shuttling well can be 
calculated in a similar way. We obtain in the end the switching rates: 



Amt fB 



^tun 

A, 



T out = D\ dBe— dAe~^ 

(7.143) 



r^shut V(B) n 
\ dBen I 

Ja„ JB 



r in = D[ I dBe~^ I dAe-^ 1 



'A; 

We recall here the equation for the current and the Fano factor for a 
dichotomous process inserting now the particular currents and switching rates 
characteristic of the SDQS coexistence regime: 



jstat _ -^shFout + -^tunFin 
Tin ~\~ r ou t 

F _£(0)_ 2 (/ sh 

-^tun) r m r ou t 

I -^shTout ~\~ ^tunTin (F m -|- F ou t) 

They represent, together with the stationary distribution 



W stat (A) = 




(7.145) 



the starting point for a quantitative comparison between the simplified model 
and the full description. 



(7.144) 
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7.3.4 Comparison 

The classical (but "noisy") limit 

The comparison between the dichotomous process model and the full descrip- 
tion for the coexistence regime is based, as for the other regimes, on the three 
investigation tools: phase-space distribution, current and current-noise. 

The phase space distribution is the most sensitive method to compare 
the model and the full description. One of the basic procedures adopted in 
the derivation of the Kramers equation (7.110) is the expansion to second 
order in the small parameters (7.96). In order to test the reliability of the 
model we simplify as much as possible the description reducing the model 
to a classical description: namely taking the zero limit for the parameter 
(sn. We realize physically this condition assuming a large temperature 

and a tunneling length A of the order of the thermal length A t h = \J 
We partially discussed this limit in the previous chapters. Also the full 
description is slightly changed, but not qualitatively: the three regimes are 
still clearly present with their characteristics. The numerical calculation is 
though based on a totally different approach 26 . 




Figure 7.15: Stationary amplitude probability distribution W for the SDQS in 
the coexistence regime. We compare the results obtained from the simplified model 
(full line) and from the full description (circles, triangles and squares). These 
results are obtained in the classical high temperature regime ksT ^> hu>. The 

amplitude is measured in units of Ath = \J 7^72" ■ ^ e mechanical damping 7 in 
units of the mechanical frequency u. The other parameter values are d = 0.05A t h 
and T = O.Olbto, A = 2A t h- 



26 The continued fraction method was applied. See for example [50] for the application 
of this method to the solution of Fokker-Planck equations. The code was developed by 
Dr. T. Novotny. 
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In figures 7.15, 7.16 and 7.17 we present respectively the results for the 
stationary Wigner function, the current and the Fano factor in the semiclas- 
sical approximation and full description. 
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Figure 7.16: Current in the coexistence regime of SDQS. Comparison between 
semianalytical and full numerical description. For the parameter values see 
Fig. 7.15. 
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Figure 7.17: Fano factor in the coexistence regime of SDQS. Comparison be- 
tween semianalytical and full numerical description. For the parameter values see 
Fig. 7.15. 

The quantum limit 

The coexistence regime in the parameters corner typically presented in this 
thesis (e.g. Fig. 5.5) is not captured by the model we introduced. The concept 
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semia nalyliGal 




• numerical 
semia nalytical 




153 



CHAPTER 7. SIMPLIFIED MODELS 



of elimination of the fast dynamics is still correct, but the small parameters 
expansion fails. The effective potential calculated from a second order ex- 
pansion still gives the position of the ring structure with reasonable accuracy 
but the overall stationary Wigner function is not reproduced due to a non 
convergent diffusion function D(A). It is clear that we need to consider 
higher order terms in the parameter (xq/X) 2 . This represents nevertheless a 
fundamental problem since it would produce terms with higher order deriva- 
tives with respect to the amplitude A in the Fokker-Planck equation and 
consequently, to our knowledge, the break down of the escape time theory. 

It has nevertheless been demonstrated with the help of the higher cumu- 
lants of the current that the dichotomous process description of the coexis- 
tence regime is valid also more in the quantum regime (A = 1.5xq), the only 
necessary condition being a separation of the ring and dot structures in the 
stationary Wigner function distribution [51]. 

Motivated by the fact that the effective potential was able to give the 
correct position of the shuttling ring we used the diffusion constant as a fitting 
parameter. The current and the Fano factor are very nicely reproduced in 
terms of a fitted diffusion constant approximately doubled with respect to 
the one calculated at zero temperature with no f 2 ^) corrections 27 (Figs. 
7.18 and 7.19). 




Figure 7.18: Current in the coexistence regime of SDQS. Comparison between 
semianalytical description with an effective fitted diffusion constant (see the text for 
explanation) and full numerical description. The parameter values are T = O.OliiJ, 
A = 2xq, d = 0.5xq, T = 0. 



Even if at this level it could appear premature we wonder if this correction can be 
attributed to an effective temperature felt by the oscillator in contact with the electrical 
system lead-dot-lead far from equilibrium. A. Armour et al. support this hypothesis in 
similar devices [40] even if in slightly different regimes. 
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Figure 7.19: Fano factor in the coexistence regime of SDQS. Comparison between 
semianalytical description with an effective fitted diffusion constant (see the text for 
explanation) and full numerical description. The parameter values are V = O.Olu;, 
A = 2xq, d = 0.5xo, T = 0. 



The full numerical description presented in this thesis is relatively fragile 
with respect to the modification of the tunneling length 28 at least for two 
reasons: the enlargement of the shuttling amplitude enlarges the number 
of oscillator states relevant for the dynamics to untractable numbers; the 
preconditioning loses effectiveness and the problem is not convergent. For 
this reason we are (still) lacking a numerical benchmark for the intermediate 
range of values of the tunneling length A > Xq but which is still non-classical 
where a second order expansion in the small parameter (^) is probably a 
good approximation. 



'We refer to changes of the tunneling length of one order of magnitude. 
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Chapter 8 
Conclusions 



We summarize in this chapter the main results presented in the thesis. The 
richness of the shuttle model still leaves many points open in the description 
and understanding of the device dynamics. We close the chapter listing some 
of the open questions that could encourage a continuation of the present work. 

8.1 Summary 

We analyzed in this thesis the dynamics of two different shuttle devices: 
the single- and the triple-dot quantum shuttles. Shuttle devices are nano- 
electromechanical systems (NEMS) in which an oscillating quantum dot can 
transfer electrons one per cycle from source to drain lead. The dynamics of 
such devices is the result of a strong interplay between their mechanical and 
electrical degrees of freedom. 

The electrical dynamics is permeated of quantum effects: electrons can 
go from the source to the quantum dot and finally to the drain lead only via 
tunneling events; the small capacitance of the oscillating dot and the corre- 
sponding large charging energy reveals the discrete nature of the electrical 
charge; finally the presence of three coupled quantum dots is also a source 
of coherent electrical dynamics in the TDQS. The small size of the device 
makes quantum effects important also for the mechanical degree of freedom. 
The resonance frequency in such small (and relatively stiff) devices can be 
in the order of GHz. Consequently the quantum of energy of the harmonic 
oscillator is comparable with other energies in the device. 

For these reasons for the SDQS we developed quantum models: we ex- 
tended for the SDQS the classical model proposed by Gorelik et al. [1] and 
we adopted for the TDQS the one already existing invented by A. Armour 
and A. MacKinnon [7]. 
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Among the different approaches for the description of the dynamics of a 
quantum system we chose the Generalized Master Equation (GME). GME's 
are particularly well-suited to study small open quantum systems (i.e. sys- 
tems with a small number of degrees of freedom in contact with resevoirs but 
still conserving some coherent dynamics) since they neglect the many (irrel- 
evant) degrees of freedom of the resevoirs still keeping track of the residual 
quantum coherencies of the small system. 

The shuttle devices are small quantum systems in contact with two kinds 
of resevoirs: the leads and the thermal bath. The strength of the coupling is 
very different and calls for a different treatment. The weak coupling between 
the thermal bath and the mechanical degree of freedom of the device justi- 
fies the perturbation theory and the quantum-optical derivation of a GME 
in Born-Markov approximation. The tunneling coupling of the shuttling de- 
vices to their electrical baths (the leads) is not weak. It sets, on the contrary, 
the time scale of the electrical dynamics that in the shuttling regime is com- 
parable with the period of the mechanical oscillations in the system. This 
strong coupling condition is realized by a tunneling amplitude modulated by 
the displacement of the quantum dot from the equilibrium position in the 
SDQS while is directly set by a model parameter in the TDQS. In order to 
handle the strong coupling we adopted the Gurvitz method for the derivation 
of the GME and extended the original formulation to take into account the 
mechanical degree of freedom. 

The shuttle dynamics has, especially in the single dot device, an appeal- 
ing simple classical interpretation and one can say that the name itself of 
"shuttle" suggests the idea of sequential and periodical loading, mechani- 
cal transport and unloading of electrons between a source and a drain lead. 
For this reason, while preserving the complete quantum treatment that we 
achieved with the GME, we wanted to keep as much as possible the intuitive 
classical picture and the possibility to handle the quantum-classical corre- 
spondence. The Wigner function distribution seemed to us a good answer 
to all these requirements. It allows a clear visualization of the numerical 
results obtained within the framework of the GME and it shows in its equa- 
tion of motion (the Klein-Kramers equation) an explicit quantum-classical 
correspondence (expansion in powers of h) [34]. 

Having the methods for the analysis of the shuttle device dynamics we 
looked for investigation tools. We first realized that a clear finger-print of 
the shuttling regime that distinguishes it from other kinds of phonon as- 
sisted tunneling is the charge-position(momentum) correlation. The Wigner 
distribution function on the device phase space (i.e. also charge resolved to 
take into account the two electrical states of the quantum dot) is a perfect 
tool to visualize this property. The only problem with the Wigner distri- 
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bution functions is that it seems very difficult (if not impossible) to access 
experimentally these functions. For this reason, with the guidance of the 
distribution function description, we investigated the electronic transport 
properties in the shuttle devices (more realistically measurable): first the 
stationary current and then the current-noise 1 . A special role in the calcula- 
tion of the current-noise for the SDQS was played by a specific characteristic 
of the GME obtained in the Gurvitz approach. In this method the informa- 
tion on the state of the leads in not completely traced out and the number 
n of electrons that passed through the device after a specific time appears in 
the equation of motion for the so-called n-resolved reduced density matrix. 

With these tools we investigated the properties of the quantum shuttle 
devices. The sharp tunneling-shuttling transition found by Gorelik et al. 
for the semiclassical model turned into a smooth crossover. We also could 
recognize the onset of the shuttling regime triggered by quantum mechanical 
noise with no external electric field acting on the oscillating quantum dot. 
The shuttle regime also revealed the characteristic current quantization (in 
our set-up a saturation due to Coulomb blockade to one electron per cycle 
in the cleanest shuttling regime) and the extremely low noise typical of a 
quasi-deterministic transport regime. With all the three investigation tools 
we adopted we could recognize three operating regimes for the shuttle device: 
tunneling, shuttling and coexistence regime. They represent for the SDQS 
the whole scenario of possible dynamics. Part of the complexity of the TDQS 
is also captured in this scheme 2 . 

The specific separation of time scales in the different regimes allowed 
us to identify the relevant variables and describe each regime by a specific 
simplified model. In the tunneling (high damping) regime the mechanical de- 
gree of freedom is almost frozen and all the features revealed by the Wigner 
distribution, the current and the current noise can be reproduced with a res- 
onant tunneling model with tunneling rates renormalized due to the movable 
quantum dot. Most of the features of the shuttling regime (self-sustained os- 
cillations, charge-position correlation and current quantization) are captured 
by a simple model derived as the zero-noise limit of the full description. Fi- 
nally for the coexistence regime we proposed a dynamical picture in terms 
of slow dichotomous switching between the tunneling and shuttling modes. 
This interpretation was mostly suggested by the presence in the stationary 
Wigner function distributions of both the characteristic features of the tun- 
neling and shuttling dynamics and by a corresponding gigantic peak in the 

1 There is already a work that, in the semiclassical regime, addresses the aspect of full 
counting statistics for these devices [5]. Results for the quantum case will soon appear 
[51]. 

2 For a further insight into this particular model see [7, 10, 12, 33]. 
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Fano factor. We based the derivation of the simplified model on the fast 
variables elimination from the Klein-Kramers equations for the Wigner func- 
tion 3 and a consequent derivation of an effective bistable potential for the 
amplitude of the dot oscillation (the relevant slow variable in this regime). 

8.2 Some open questions 

During the developing of the present work many questions arose about pos- 
sible phenomena to investigate in different shuttle devices or other possible 
regimes that the systems we investigate could exhibit in some other param- 
eter corner that we did not reach. This questions could be of interest for 
example in understanding or proposing future experiments. We list them 
here as a spur for possible continuation of the present work sure that some 
of the theoretical tools presented in this thesis are definitely useful for those 
tasks. 

• The mechanical bath introduces decoherence in the mechanical degree 
of freedom, the leads in the electrical. Correlations are left (in the 
SDQS) only in the off-diagonal elements of the charge resolved density 
matrices. Can we imagine an electromechanical system that instead ex- 
hibits quantum coherence in the mechanical degree of freedom? Could 
we use the electromechanical interaction to "pump" continuously co- 
herence in the device? 

• The stationary solution for the GME is incoherent. Is there any possi- 
bility to reduce the GME to an effective Master equation? We imagine 
that especially in the coexistence regime this could help to avoid the 
Wigner function small parameter expansion for the calculation of the 
effective bistable potential and lead directly to the calculation of the 
rates for the slow switching process. This would make the calculation 
of the parameters for the dichotomous process description much more 
general. 

• Taking into account the spin degrees of freedom and the Coulomb block- 
ade the incoming bare rate should be considered as double with respect 
to the out going rate. This fact does not change qualitatively the pic- 
ture since the switching is, in first approximation a non-dissipative 
process and the electrostatic force is a small perturbation of the oscil- 
lator restoring force. The mechanical trajectories are not significantly 

3 The analytical derivation of the effective potential is an extension of the work done 
by Fedorets et al. [9]. 
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changed (maybe slightly decreased in radius) and only the symme- 
try between the two charged resolved Wigner functions is modified in 
favour of W\\. the dot is charged for a longer time. But what happens 
if we consider a spin dynamics on the dot due for example to a coupling 
to an external magnetic field ? 

• At high bare injection rates (r ^> a;) the shuttle device is more stable 
and hardly leaves the tunneling regime. The mechanical motion can- 
not in fact follow the fast electrical forcing and the effect is generally 
uncorrelated resonant tunneling. Nevertheless at extremely low damp- 
ing the shuttling regime arises. The transition is very different from 
the small injection rate shuttling instability: it is very smooth and the 
ring structure gradually emerges from the central fuzzy spot typical for 
the tunneling regime. The general question is: can we draw a phase 
diagram for shuttling devices as a function of their parameters? 

• A big issue regards the treatment of the bias in the quantum descrip- 
tion. All our equations were derived in the limit of high external bias. 
We did not have access to this parameter and we used for this reason 
the mechanical damping as a control parameter. Is it possible to derive 
a GME at low bias and for strong couplings to the leads? Probably 
the dynamics will get non-Markovian and all kinds of coherent trans- 
port effects will also enter the game. We know from the semiclassical 
treatment that the shuttling instability should be present. Can we 
understand something new from the quantum treatment? 



4 A partial answer to this questions is contained in the very recent paper by Fedorets 
et al. [47]. 
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